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I .  INTRODUCTION 


It  is  quite  important  to  classify  the  ships  according  to 
their  types.  Classification  can  be  done  in  a  number  of 
different  ways.  The  simplest  is  the  visual  method  that  is 
prone  to  error.  Other  methods  for  classification  are  the 
Fourier  Coefficient  method  and  the  B-spline  Coefficient 
method.  In  these  methods  the  necessary  information  for  clas¬ 
sification  can  be  obtained  from  the  superstructure  profile. 

The  Fourier  Coefficient  method  samples  the  superstruc¬ 
ture  profile  at  every  chosen  points.  The  function  values  at 
the  sampling  points  are  transformed  into  the  spatial  compo¬ 
nents.  The  logarithm  of  the  magnitude  of  these  components 
is  plotted  and  compared  with  the  standard  plot  to  recognize 
the  type  of  the  ship. 

In  the  B-spline  Coefficient  method,  the  spline  coeffi¬ 
cients  along  the  X  axis  and  the  Y  axis  are  used  to  recon¬ 
struct  the  superstructure  profile.  The  shape  of  the  curve  of 
the  spline  coefficients  is  in  some  way  similar  to  the  shape 
(the  position  of  the  lumps)  of  the  superstructure  profile. 
The  ship  classification  may  be  achieved  by  recognizing  the 
beginning,  the  peak,  and  the  area  of  the  lumps. 


II.  PREPROCESSING 


Preprocessing  is  the  procedure  to  obtain  the  superstruc¬ 
ture  profile  of  a  ship  from  the  IR  image.  Then,  Fourier 
coefficient  or  spline  coefficient  methods  can  be  used.  The 
details  of  the  preprocessing  procedures  are  a  follows. 

A.  DATA  COLLECTION 

The  data  consists  of  the  IR  image  of  eight  different 
types  of  ships. 

1.  DD  -  Destroyer;  "HALL"  class. 

2.  Container;  The  class  is  unidentified. 

3.  Freighter;  The  class  is  unidentified. 

4.  AOR  -  Replenishment  oiler;  "WICHITA"  class. 

5.  LST  -  Tank  landing  ships;  "NEWPORT"  class. 

6.  FF  -  Frigate;  "GARCIA"  class. 

7.  CGN  -  Guided  missile  cruiser  (Nuclear  propulsion); 
"BAINBRIDGE"  class. 

8.  DDG  -  Guided  missile  destroyer;  "CHARLES  F.  ADAMS" 
class . 

Tb«g<»  -'-'ages  are  taken  from  an  aircraft  which  is  flown 
at  a  ........  .  500  feet  with  a  speed  of  400  knots  toward  the 

side  of  the  ship.  The  aspect  angles  for  these  images  are  90 
degrees  which  may  be  slightly  off  in  some  images.  The  inac¬ 
curacy  of  the  aspect  angle  arises  from  the  fact  that  the 
photos  are  taken  while  the  aeroplane  keeps  on  moving.  All 
data  of  the  images  are  stored  in  a  digital  magnetic  tape 
with  256  by  64  bytes  per  image.  Thus,  the  number  of  bytes 
required  for  each  image  is  16384  and  each  byte  represents 
the  intensity  of  a  pixel.  For  each  record  of  the  image,  a 
label  is  coded  in  the  last  8  bytes  as  follows: 

1.  Byte  (16377)  =  Run  number  in  each  flight  which  passes 
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Figure  2.16  Contour  Image  of  a  CGN  at 
a  Range  of  45000  feet. 


H.  CLOSING  OPERATION 

Some  results  from  the  superstructure  extraction  process 
are  discontinuous  because  the  gray  level  of  those  areas  of 
the  structure  is  less  than  the  threshold  value.  If  we 
decrease  the  threshold  value,  the  details  of  the  superstruc¬ 
ture  are  effected.  It  is  necessary  to  use  the  "Closing" 
operation  which  consists  of  the  "dilation"  process  to  smooth 
the  superstructure  profile  used  the  "Erosion"  process.  All 
direction  are  dilated  similarly  which  causes  an  smoothing 
effect  on  the  edges.  The  superstructure  increases  in  total 
area.  Then,  use  the  "Erosion"  process  to  shrink  (subtract) 
the  dilated  part  in  all  directions,  thus  obtain  the  smoothed 
superstructure  with  the  appropriate  size.  The  Closing 
process  employs  the  dilation  process  which  yields  an  output 
image  from  an  input  image.  The  Dilation  results  is  shown  in 
Figure  2.17. 


image  is  shown  in  Figure  2.16.  If  the  superstructure  in 
Figure  2.9  has  wide  edge,  we  can  not  find  the  contour  image. 
We  have  to  use  the  additional  step  which  is  called  the 
"Closing"  operation. 


Figure  2.10  The  3  by  3  Kernel. 


Starting  from  the  leftmost  point  in  the  superstructure 
image  in  Figure  2.9  the  contour  profile  tracing  is  accom¬ 
plished  by  examining  the  neighbors  of  a  3  by  3  kernel 
located  at  the  curser  position  as  shown  in  Figure  2.10.  The 
curser  is  moved  along  the  profile.  All  successive  positions 
of  the  curser  constitute  the  contour  profile  of  the  ship. 
The  testing  procedure  is  explained  below. 

1.  Initialize  the  curser  position  to  the  beginning  of 
the  thresholded  image  with  the  gray  level  of  255  and 
the  estimate  direction  of  the  next  position. 

2.  Check  for  reaching  the  end  of  the  profile,  if  it  is 
at  the  ending  of  the  profile  then  stop,  if  not  go 
to  3 . 

3.  Check  for  the  estimate  to  see  whether  it  is  in  the 
direction  of.  North,  East,  South,  or  West.  If  the 
direction  it  is  North  then  go  to  4,  if  it  is  East 
then  go  to  5,  if  it  is  South  then  go  to  6,  and  if  it 
is  West  then  go  to  7. 

4.  Determine  the  next  position  with  the  gray  value  of 
255  in  the  counter-clockwise  direction  as  shown  in 
Figure  2.11.  Store  the  position  found  and  move  the 
curser  to  that  position.  Update  the  estimate 
direction  to  the  one  last  found;  then  go  to  2. 

5.  The  procedure  is  the  same  as  that  in  step  4  except 
search  pattern  is  shown  in  Figure  2.12. 

6.  The  procedure  is  the  same  as  that  in  step  4  except 
search  pattern  is  shown  in  Figure  2.13. 

7.  The  procedure  is  the  same  as  that  in  step  4  except 
search  pattern  is  shown  in  Figure  2.14. 

The  flow  chart  of  the  procedure  is  shown  in  Figure 
2.15,  and  the  detail  of  each  procedure  are  included  in 
Appendix  B.  The  testing  procedure  have  to  be  performed  in 
such  a  maner  that  the  resulting  contour  is  a  good  represen¬ 
tation  of  the  superstructure  line.  The  result  of  the  contour 
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slope  line  to  zero.  Hence,  it  results  in  a  superstructure  of 
the  ship  as  shown  in  Figure  2.9.  However,  in  some  images, 
there  is  a  lot  of  noise  in  the  background  which  cause  diffi¬ 
culty  in  locating  the  bow  and  stern.  Under  these  circum¬ 
stances,  the  dimensions  of  a  ship  are  estimated  by  trial  and 
error  method.  Then  the  gray  level  which  lies  outside  is  set 
to  zero  and  an  estimate  of  the  bow  and  stern  slope  can  be 
made. 


Figure  2.9  Superstructure  Profile  of  Figure  2.8. 


G.  CONTOUR  FOLLOWING 

The  noise  in  the  background  of  the  original  image, 
yields  wide  edge  structure.  In  considering  this  factor,  the 
contour  following  procedure  tracking  the  inner  part  of  the 
image  in  Figure  2.9  gives  the  superstructure  profile.  The 
objective  of  this  contour  following  is  to  describe  the  bow 
and  stern  points  of  the  ship,  the  direction,  and  the  posi¬ 
tion  of  the  edge  of  the  superstructure.  The  contour  tracing 
is  done  in  the  counter-clockwise  direction  which  compares 
pixel  value  of  0  or  255  in  a  3  by  3  matrix  in  the  following 


manner. 


E.  HOW  TO  EXTRACT  THE  PROFILE 


The  edge  image  of  the  guided  missile  cruiser  shows 
little  variations  in  the  gray  level  as  shown  in  Figure  2.3. 
These  variations  are  caused  by  the  noise  in  the  original 
image.  In  this  case,  the  choice  of  the  threshold  value  is 
based  upon  both  the  histogram  and  the  cumulative  distribu¬ 
tion  of  the  edge  image  so  that  it  contains  90  %  of  the 
pixels.  Therefor,  the  chosen  gray  level  is  110  and  the 
result  is  shown  in  Figure  2.8. 


Figure  2.8  Silhouette  of  a  CGN  in  Figure  2.3. 


F.  HOW  TO  OBTAIN  THE  SUPERSTRUCTURE 

The  original  image  .  of  the  ship  is  taken  from  the  aero¬ 
plane  with  different  displacement  from  waterline  to  the 
superstructure  in  a  rough  sea,  so  that  the  profile  of  the 
ship  with  respect  to  the  sea  surface  varies.  Thus,  we  have 
to  eliminate  some  information  in  the  image  of  the  ship  by 
considering  the  superstructure  only.  In  considering  the 
overall  ship  structure,  it  is  obvious  that  the  largest 
distance  is  between  the  bow  and  stern  span.  Therefore,  the 
bow  and  stern  points  are  located  first  in  the  program  as 
shown  in  Appendix  A.  Then  we  consider  the  slope  of  the  bow 
and  stern  of  the  ship  and  set  all  the  gray  values  below  the 


D.  EDGE  THRESHOLD  STRATEGIES 

Edge  threshold  strategies  are  used  to  extract  the  edge 
profiles  from  the  Sobel  results.  In  this  case,  we  use  only 
the  pronounced  value  of  an  edge  element  atxif  g(x)  is 
greater  than  certain  threshold  value  [Ref.  2). 


=  gfoy) 


G(x,y)  < 


if  g(x  y)^.  threshold 


othe  rwi se 


(2.5) 


To  increase  the  contrast  of  the  image  to  a  silhouette 
form,  G(x,y)  is  defined  as 


G(x,y) 


f  =255 


1=0 


if  9(x;y) ^threshold 


otherwise 


(2.6) 


The  choice  of  the  threshold  value  is  based  upon  the 
histogram  of  the  edge  image  as  shown  in  Figure  2.6.  The 
estimated  critical  gray  level  is  chosen  so  that  a  majority 
number  of  the  pixels  with  value  between  0  to  255  will  fall 
below  the  critical  value.  Alternatively,  histogram  equaliza¬ 
tion  may  be  used  to  determine  the  desired  threshold  level  as 
shown  in  Figure  2.7.  In  this  case,  trial  and  error  method 
was  used  in  conjunction  with  the  above  method  to  obtain  the 
threshold  so  that  the  correct  profile  is  ascertained. 
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Figure  2.1  Sobel  Operator. 


C.  THE  USE  OF  A  SOBEL  OPERATOR 

When  a  Sobel  operator  is  used  at  the  edge  of  the  image 
frame,  the  pixel  level  which  lies  out  of  the  frame  will  be 
set  equal  to  that  of  the  adjacent  pixel  within  the  frame. 
The  original  image  of  a  guided  missile  cruiser  is  shown  in 
Figure  2.2.  Since  the  result  of  the  Sobel  operator  is  numer¬ 
ically  greater  than  8  bits  range,  we  have  to  rescale  the 
result  back  to  8  bits  range.  This  is  achieved  by  deter¬ 
mining  the  maximum  and  the  minimum  of  the  gray  level.  They 
are  then  used  to  rescale  the  gray  level  in  the  results  of 
the  Sobel  operator  as  shown  in  Figure  2.3.  In  some 
instances,  the  original  image  is  very  poor  as  shown  in 
Figure  2.4.  Attempts  to  determine  the  edge  of  this  image 
failed  as  shown  in  Figure  2.5. 


The  basis  vector  for  all  directions  are  (a-2b+c), 
(g-2h+i),  (a-2d+g),  and  (c-2f+i).  Furthermore,  each  of  the 
basis  vector  is  convolved  with  the  image  as  follows: 
along  the  x  axis 


d*=  [f(x-  l,y- 1)+ 2f(x,y-i)+f(x+l,y-i)] ' 


~[f(x~U  y-H)+ 2f(x;y+ 1) +f(x+i,y+ !)] 


(2.2) 


along  the  y  axis 


dy  =  [f(x+l,y-|)+2f(x  +  U)+f(x+li  y+ 1 )] 

“  [f  (x  -  l,y-  0+  2f  (x- y ) + f  (x-  |,y+ 1 )] 


(2.3) 


Since  the  magnitude  of  the  resulting  vector  is  the  abso¬ 
lute  value  of  the  convolved  results,  the  edge  magnitude 
S(x, y )  [Ref.  1], 


S(x,y)  =  (d/+d j)'!2 


(2.4) 


Note  that  the  Sobel  operator  does  not  use  the  gray  level 
at  the  position  (x,y).  The  advantage  of  using  a  Sobel  oper¬ 
ator  over  others  is  that  the  resulting  edge  is  smoother  due 
to  a  3  by  3  matrix  approach.  If  we  compare  the  Sobel  oper¬ 
ator  with  the  Laplacian  operator,  it  is  seen  that  the  Sobel 
operator  using  the  four  basis  vector  as  shown  above  will 
provide  more  accurate  reading  because  of  noise  reduction  in 
the  original  image.  Hence,  The  Sobel  operator  is  often  used 
in  the  preprocessing  operation. 


*  •  *  •  *  * 
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the  ship. 

2.  Byte  (16378)  =  Video  tape  time  code  when  the  data  is 
taken;  in  minutes. 

3.  Byte  (16381)  =  Video  tape  time  code;  in  seconds. 

4.  Byte  (16382)  =  Video  tape  time  code;  in  thirtieth  of 
a  second. 

5.  Byte  (16379)  =  Range  in  kilo-feet  which  is  the 
distance  measured  from  the  radar  it  may  have  an 


from 


error  1  to  2  kilo-feet. 

6.  Byte  (16380)  =  Aspect  angle;  degrees  from  the  bow 
of  the  ship. 

7.  Byte  (16383)  =  Ship  class. 

8.  Byte  (16384)  =  ID,  1  =  for  training,  2,  3,  4,  5  =  for 
testing. 

The  run  number  and  the  time  code  together  uniquely 
define  a  specific  image  that  represents  a  single  TV  frame 
with  no  averaging.  In  addition,  the  time  interval  between 
the  end  of  one  image  to  the  beginning  of  the  other  is 
approximately  1.5  seconds.  Also,  there  are  inherent  random 
noise  in  the  record  which  arises  from  the  photo  instrument 
and  the  process  of  storing  them  on  to  the  digital  magnetic 
tape. 


B.  SOBEL  OPERATOR 

The  Sobel  Operator  technique  is  used  to  find  the  edge. 
To  determine  the  edge,  the  Sobel  Operator  uses  the  differ¬ 
ence  of  gray  levels  of  the  pixels  in  a  3  by  3  matrix  as 
shown  in  Figure  2.1. 

a,b,c,d,e,f,g,h,  and  i  are  the  values  of  the  gray 
levels  at  the  position  of  (x-l,y),  (x,y-l),  (x-l,y),  (x,y), 
(x+l,y),  (x-l,y+l),  (x,y+l),  and  (x+l,y+l)  respectively.  The 
Laplacian  estimate  is  defined  as 


~  f(x,y)  - f(x+  l,  y)  -[f(x+  l,y)-f(x+  2,y)] 


(2.1) 


Figure  2.17  The  Process  of  Dilation. 


We  examine  the  gray  value  of  each  pixel  in  the  original 
image.  Begining  from  the  pixel  at  the  first  row  and  the 
first  column.  The  procedures  are  the  following. 

1.  Considering  one  pixel  in  the  original  image  with  the 
kernel  (B)  centered  there.  If  at  least  one  pixel  in 
the  kernel  has  a  value  of  255,  we  let  the  gray  level 
of  the  output  image  at  the  center  of  the  kernel  be 
255. 

2.  We  shift  the  center  of  the  kernel  1  column  to  the 
right.  Then  following  the  same  procedure  as  in  step  1 
until  the  last  column  is  reached. 

3.  We  shift  the  position  of  the  kernel  to  the  next  row 
and  starting  from  the  first  column.  Then,  following 
the  same  procedure  as  in  step  1  and  step  2  until  the 
last  row  and  the  last  column  is  reached. 

The  result  obtained  from  the  dilation  process  is  an 
image  with  enlarged  structure.  The  second  procedure  in  the 
Closing  operation  is  the  Erosion  process.  The  Erosion 
process  perform  the  same  procedures  as  the  dilation  process 
except  for  step  1.  If  every  pixel  in  the  kernel  are  255,  we 
let  the  gray  level  in  the  new  image  at  the  center  of  the 
kernel  be  255.  Otherwise,  it  will  be  0.  The  output  image 
obtained  will  have  smooth  edge  with  minor  change  occurring 
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in  the  edge  detail  as  shown  in  Figure  2.18.  In  this  case,  we 
use  an  kernel  of  size  3  by  3 .  If  we  increase  the  size  of  the 
kernel,  the  details  of  the  image  are  decreased. 


Figure  2.18  The  Profile  after  Dilation  and  Erosion 

(Closing  Process). 


I.  PROFILE  ROTATION 

Often  in  the  contour  image,  the  first  and  the  last  point 
of  the  superstructure  are  not  at  the  same  horizontal  level. 
We  have  to  rotate  the  contour  image  by  setting  the  two 
points  to  the  same  horizontal  level.  How  to  rotate  it  from 
the  Y,  X  axis  to  the  Y' ,  X'  axis  is  shown  in  Figure  2.19. 

If  the  angle  value  Q  is  positive,  it  will  be 
fx'l  _  [cos Q  -sing!  [X] 

[y'J  [sin0  COS0J  LYJ 

If  the  angle  value  Q  is  negative,  it  will  be 
[X'l  =  [  cos Q  sing!  [Xl 
[y'J  It  sing  cosgj  [yJ 

For  some  of  the  contour  profile,  part  of  the  profile 
after  rotation  will  be  out  of  the  image  frame.  Then  we  have 
to  shift  this  contour  down  by  20  pixels  position.  The 
rotated  profile  is  shown  in  Figure  2.20. 
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Figure  2.19  Rotation  Process. 


III.  FOURIER  COEFFICIENT  METHOD 

We  have  obtained  the  ship  profile  in  the  previous 
chapter.  To  extract  features  out  of  this  profile  for  clas¬ 
sification  purposes,  we  will  use  the  Fourier  Transform 
method . 

The  Fourier  transform  of  the  ship  profile  showed  that 
the  transform  coefficients  depend  upon  the  ship's  dimen¬ 
sions,  its  superstructure,  and  the  distance  between  the 
camera  and  the  ship.  If  the  profile  f(x)  is  a  discrete 
function  with  128  sample  points.  The  discrete  Fourier 
transform  can  be  written  as 

F(u)  =l^f(x)exp(-j27Tu  x/N)  (3  1} 

For  u,  x  =  0,  1,  2,  . .  N-l.  N  is  the  total  number  of 

samples. 

If  the  direct  calculation  of  the  discrete  Fourier  trans¬ 
form  is  chosen,  the  number  of  complex  multiplication  ar_d 
addition  will  be  equal  to  N  ;  i.e.  to  obtain  F(0)  would 
require  complex  multiplication  and  addition  N  times.  In 
order  to  reduce  the  computation,  we  use  the  fast  Fourier 
transform  algorithm.  Thus,  equation  3.1  [Ref.  2]  can  be 
separated  into  Few|  (u)  and  F^  (u) 


W  =  e  x  p(-j  27T/M) 

M  =Ji 
2 

F(U>=l[W“>  +  F.<>)W2*] 


F(u+M)=1[f„,»-F0<(»^J 

Using  this  method,  we  have  reduced  the  total  number  of 
complex  multiplication  and  addition  to  Nlog  N.  In  this  case, 
we  have  N  =  128,  thus,  the  total  number  of  complex  multipli¬ 
cation  and  addition  will  be  896. 

First,  we  divide  the  rotated  profile  image  into  128 
divisions.  Since  the  distance  of  each  division  is  equal  to 
the  amount  of  the  pixel  between  the  bow  and  the  stern  of  the 
ship  divided  by  128,  we  use  the  distance  perpendicular  from 
the  horizontal  line  between  bow  and  stern  to  the  highest 
point  of  the  superstructure  as  the  sampled  values.  The 
result  of  the  fast  Fourier  transform  is  a  complex  number. 
Then  we  use 


(3.4) 

(3.5) 

(3.6) 

(3.7) 


M(u)  =  lo  g[|+G(u)]  0.8) 

G(u)  is  the  magnitude  of  F(u).  A  value  of  1  is  add  to  the 
magnitude  to  avoid  negative  logarithm  result,  the  results 
obtained  are  as  shown  in  Table  I  and 

1.  Figure  3.1  -  Destroyer 

2.  Figure  3.2  -  Container 

3.  Figure  3.3  -  Freighter 

4.  Figure  3.4  -  Replenishment  oiler 

5.  Figure  3.5  -  Tank  landing  ship 
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6.  Figure  3.6  -  Frigate 

7.  Figure  3.7  -  Guided  missile  cruiser 

8.  Figure  3.8  -  Guided  missile  destroyer 

On  the  results  minor  difference  of  the  shape  of  the 
Fourier  coefficients  can  be  noticed  visually.  But  it  is 
difficult  to  implement  a  program  to  detect  these  minor 
difference  in  shape.  Therefore,  a  Second  method  was  used  to 
handle  this  problem. 


Figure  3.1  Logarithmic  Magnitude  of  the  Fourier  Transform  of  a  DD. 


Magnitude  of  the  Fourie 


Figure  3.3  Logarithmic  Magnitude  of  the  Fourier  Transform  of  a  Freighter 


Logarithmic  Magnitude  of  the  Fourier  Transform  of  a  AOR. 


Figure  3.5  Logarithmic  Magnitude  of  the  Fourier  Transform  of  a  LST 
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Plots  of  the  reconstructed  profile  of  three  ships  of 
the  CGN  class  using  different  value  of  S,  result  in 
different  number  of  knot  positions  and  are  shown  for  compar¬ 
ison  in  Figure  4.3,  4.4  and  Figure  4.5.  The  original 
samples  are  plotted  in  solid  curve,  the  reconstructed  curve 
is  plotted  in  dashed  line  in  Figure  4.3. 
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Limitation  on  B-spline  Coefficient  Determination 


If  the  value  of  NEST  is  too  small,  the  user  will 
receive  error  code  IER  and  the  number  of  knot  positions  will 
appear  to  be  dense  in  the  first  part  of  the  curve.  While 
sparse  in  the  other  part.  The  example  of  an  incorrect  selec¬ 
tion  of  NEST  is  shown  in  Figure  4.2. 


Figure  4.2  Knot  Selection  if  the  NEST  Parameter  is  too  Small 

Another  problem  which  causes  difficulty  in  program¬ 
ming  is  that  the  main  program  is  in  PASCAL  while  the  subrou¬ 
tine  is  in  FORTRAN.  As  for  the  PARAM  program  in  FORTRAN,  the 
array  index  value  starts  from  1,  while  ,  for  the  user  PASCAL 
program  it  starts  from  0.  Therefore,  in  linking  the  main 
program  to  the  subroutine,  one  has  to  keep  in  mind  of  the 
difference.  In  addition,  the  FORTRAN  programming  logic 
structure  is  so  complicated  that,  when  an  error  occurs,  it 
is  very  difficult  to  debug  and  locate  the  error. 

The  values  of  Cx  and  Cy  depend  upon  uneven  knot 
positions,  and  they  contribute  controlling  effect  to  the 
reconstructed  curve  which  would  be  both  smooth  and  close  to 
the  original  curve.  In  running  the  B-spline  approximation 
program  for  the  first  time,  the  values  of  Cx  and  Cy  of  the 
last  4  knots  at  the  right  end  are  zeros.  However,  in  running 
it  again,  increase  the  value  of  S  did  not  yield  zero  values 
of  Cx  and  Cy  at  those  points.  This,  nevertheless,  has  no 
effect  on  the  reconstruction  of  the  curve. 
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Figure  4.1  Plot  of  the  Original  Data  and  the  Approximate 

with  Free  Knot. 

C.  TO  DETERMINE  THE  KNOT  POSITION  AND  THE  B-SPLINE 
COEFFICIENTS 

The  approximate  smoothing  factor  in  eq(4.7)  is  to  be 
calculated  before  using  the  subroutine  "PARAM"  [Appendix  Cj. 
In  the  subroutine,  there  is  a  check  on  S  factor  every  time 
it  is  run.  If  the  S  input  values  do  not  satisfy  or  meet  the 
criteria,  the  program  will  return  with  a  error  code  IER . 

There  is  a  parameter  NEST  which  is  set  to  a  constant 
value.  This  value  determine  the  dimension  of  the  knot  posi¬ 
tions  which  relates  to  the  array  T(NEST),  Cx(l..NEST)  and 
Cy(l..NEST).  The  NEST  is  an  overestimate  of  the  dimension  of 
the  arrays  set  by  the  user.  The  limitation  on  the  value  of 
NEST  for  subroutine  "PARAM"  are  as  follows  [Appendix  C] 

1.  2k+2  NEST  M+k+1 

2.  Typically,  value  of  NEST  M/2 

where  M  is  the  number  of  the  total  sampling  points  and  k=3 
is  the  highest  order  in  the  B-spline  function. 

The  subroutine  "PARAM"  produces  several  outputs  as,  N 
the  total  number  of  knot  positions,  T(N)  the  array  of  the 
value  of  knot  positions  in  Z  parameter,  Cx(N)  and  Cy(N), 
B-spline  coefficient  of  X  functions  and  Y  functions  for  knot 
positions  defined  in  array  T(N). 


W  can  be  equal  to  one  and  the  S  can  be  determined  by  the 
trial  and  error  method. 


1.  Interpolation  and  Approximation  Using  B-spline 

Function 

There  is  a  difference  between  the  interpolation  and 
the  approximation.  In  interpolation,  the  number  of  knot 
positions  required  are  equal  to  the  number  of  sampling 
points  and  the  value  of  S  in  eq(4.7)  is  small.  In  this  case, 
the  computation  time  required  will  increase  tremendously. 
Whereas,  in  approximation,  it  is  not  of  great  concern  that 
the  approximated  value  at  a  position  be  the  same  value  of 
the  sampling  point,  and  most  of  the  information  still  remain 
in  the  original  curve.  Thus,  decreasing  the  value  of  S  will 
result  in  the  large  number  of  the  knot  positions  and  the 
final  curve  obtained  will  be  similar  to  the  original  curve. 

The  justification  for  using  approximation  rather 
than  interpolation  approach  is  that  though  the  resulting 
curve  may  not  be  the  best  fitting  curve,  it  is  smooth  and 
close  enough  to  the  original.  The  approximation  spline  func¬ 
tion  has  less  number  of  knots  than  the  number  of  samples, 
which  reduces  the  total  processing  time.  In  approximation 
approach,  when  the  appropriate  choice  of  S  value  is  made,  it 
will,  in  turn,  generate  sufficient  number  of  knot  positions 
required  to  provide  a  close  approximation  to  the  original 
curve.  The  ratio  of  sampling  point  to  spline  coefficient  is 
10  to  1  as  shown  in  Figure  4.1.  In  Figure  4.1,  a  guided 
missile  cruiser  ship  at  a  distance  of  45000  feet,  with  290 
sampling  points  and  the  associated  B-spline  knots  are  shown. 
The  original  samples  are  plotted  in  solid  curve.  The  knot 
position  are  show  by  small  circles  in  Figure  4.1.  After 
B-spline  approximation,  the  number  of  the  knots  are  reduced 
to  33.  Hence,  this  technique  is  essentially  a  kind  of  data 
compression  technique. 
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lit: 


Alf(Z;,Z ,>/  ,  .  .  .  ,Zf-+jf  )f (t)  stands  for  the  1-th  divided 
difference  of  the  function  f(t)  on  the  point 
where  t  is  the  position  values  of  knots  in  term  of  Z  param¬ 
eter.  Z  parameter  is  defined  as  follows 


2(1)  =2(1- 1)+  [(X(0-X(I-ljf+(Y(l)-Y(I-l))2/2 


(4.5) 


2(0)  =  0 


(4.6) 


where  I  is  the  number  of  the  sampling  point  and 
X  — *1 ,2,3,  .  .  .  ,  m . 

Second,  the  smoothing  is  subjected  to  a  constraint 


6(c)  < 


(4.7) 


where  S  is  the  smoothing  factor,  6(C)  is  the  weighted  sum  of 
the  square  residuals  defined  as 


6(e)  =Iwy[ry-jc,-M.<+/(x.)]‘ 


(4.8) 


Xy  , Yy  are  the  values  of  X  and  Y  at  Z  parameter  of  the 
sampling  points;  Wy  is  a  weighting  factor  for  all  sampling 
points  (Ref.  3]  defined  as 


W;  =  (Sy/ 


(4.9) 


where  Y^-  is  an  estimate  of  the  standard  deviation  of  Y y  . 
Then,  the  value  of  S  is  in  the  range  m+V2m.  if  nothing  is 
known  about  the  statistical  standard  deviation  of  Yy  ,  each 


B.  B- SPLINE  APPROXIMATION  WITH  FREE  KNOT 


As  mention  earlier,  the  B-spline  function  is  used  in  the 
spline  approximation  with  free  knot.  The  free  knots  is  some¬ 
times  called  uneven  or  irregular  knots;  that  is  no  fixed 
number  of  knots  are  used.  Furthermore,  the  spline  position 
need  not  be  on  the  original  curve  so  as  to  minimize  the 
number  of  knot  position  while  preserving  most  of  the  detail's 
of  the  original  curve. 

First,  minimize  the  value  of  the  lack  of  the  smoothness 
7](c)  defined  as  [Ref.  3] 

7?(C>  =  Z  (  E  o,..C;  ) 

jz!  >--k  (4.  !) 


where  Ci  is  a  coefficient  of  spline  at  the  knot  position, 
a ij  is  defined  as  follows 

0,7  =  m!‘UV+0)-^+/(v-0)  (4.2) 

where  is  the  normalize  B-spline  function  and  defined 

as 

^ /';*+/ (x)  “  i\‘+  (t/ / •  •  •  )\-+  1 ) 9* (t i  x ) 

(4.3) 

and  G*(t;x),  t  are  defined  as  follows 

f=(t-x)+=  (t-x)*  if  »>x 


if  t<  x 


(4.4) 


IV.  B-SPLINE  COEFFICIENT  METHOD 


Using  B-spline  coefficient  is  another  method  to  describe 
a  ship  profile.  This  method  uses  the  B-spline  coefficients 
to  determine  the  beginning,  peak,  and  area  of  lumps  which 
contain  significant  information  about  the  type  of  the  ships. 
The  comparisons  of  the  knot  position  (in  parametric  value) 
from  the  midships  to  the  peak  or  beginning  of  the  lump,  can 
be  very  helpful.  Different  ships  will  have  different  lump 
positions  and  sizes. 

A.  BACKGROUND 

A  spline  function  is  a  piecewise  polynomial  used  to 
interpolate  points.  This  kind  of  curve  is  smooth  and  the 
discontinuities  in  its  k-th  derivative  is  as  small  as 
possible.  In  this  case,  Cubic  spline  was  used,  where  the 
first  and  second  derivatives  for  any  set  of  interpolating 
points  are  continuous,  while  the  third  derivative  may  be 
discontinuous.  The  reason  for  using  Cubic  spline  is  to  keep 
the  third  derivative  discontinuity  as  small  as  possible  and 
the  curve  as  smooth  as  possible.  The  B-spline  calculation 
procedure  is  also  very  stable.  In  our  case  the  order  of 
spline  used  is  3,  while  the  1-st  and  2-nd  derivative  are 
continuous.  The  choice  of  4-th  order  spline  function  is  due 
to  the  fact  that  it  is  generally  sufficient  for  most  ship 
profile  curves.  Discontinuity,  in  a  sense,  may  be  stated  as 
the  jumps  of  the  third  derivative,  which  is  the  means  to 
control  smoothness  of  tl  a  connecting  pieces. 


of  the  field  of  view  is  accurate,  we  can  estimate  the  number 
of  pixel  in. an  object  image.  Then  we  can  classify  the  object 
by  comparing  the  number  of  pixel  of  the  original  image  with 
that  of  the  test  image.  Some  known  system  parameters  can 
help  to  determine  the  range  of  the  target  ship.  The  problem 
is  that  existing  errors  in  the  system  parameter  causes  in 
the  larger  errors  in  the  estimated  range.  Consequently,  they 
are  not  very  useful  in  classifying  the  ships. 


TABLE  II 

Range  Estimation 


CLASS 

I (pixel) 

MEASURE 

1  0  (ft) 

|  KNOWN 

|  R  (kft)  | 

| CALCULATE] 

R(kft)  |  (R-R)IOO 

1  . 

RADAR  0IS|  R 

1 

D01 

96 

!  418 

|  21.766  I 

77 

1 

|  71.71 

002 

80 

|  418 

I  26.119  i 

85 

|  69.27 

A0R1 

107 

I  659 

I  30.787  | 

78 

|  60.53 

A0R2 

96 

i  659 

|  34.315  | 

85 

j  59.63 

LST1 

176 

|  522.3 

|  14.834  I 

51 

|  70.91 

LST2 

134 

|  522.3 

I  19.484  | 

57 

|  65.82 

j 

CCN1 

147 

|  565 

I  19.213  | 

45 

|  57.31 

CGN2 

119 

|  565 

|  23.734  I 

55 

|  56.85 

DDC1 

126 

1  437 

I  17.337  | 

47 

3  63.11 

DDG2 

90 

|  437 

I  24.247  | 

64 

|  62.08 

DD1.DD2  *  Destroyer  *t  range  79QOO  and  83000  feet. 

A0R1 , A0R2  =  Replenishment  oiler  at  range  78000  and  88000  feet. 
LST1.LST2  *  Tank  landing  ship  at  range  51000  and  62000  feet. 
CCN1.CCN2  *  Guided  missile  cruiser  at  range  45000  and  64000  feet. 
0DC1.DDG2  *  Guided  missile  destroyer  at  range  41000  and  64000  feet. 


The  error  distance  between  the  estimated  distance 
and  calculated  distance  in  percentage  is  ( (R  -  R)/R)  100. 

2 .  Remark 

Calculated  distance  error  in  R  may  come  either  from 
the  pixel  measurement  in  an  image  or  from  the  angular  reso¬ 
lution  estimation.  The  distance  that  is  stored  in  the  image 
label  has  an  error  of  about  1  to  2  kilo-feet.  If  the  angle 
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2  tana 
2 

Assume  that  the  angle  resolution  of  the  pixel  of  the 
field  of  view  of  the  camera  is  equal  to  0.2E-3  radian  per 
pixel.  The  number  of  pixel  of  the  frame  is  equal  to  256. 
The  size  of  an  image  is  I  pixels.  The  dimension  of  the 
object  D  in  feet  is  known.  The  field  of  view  in  angle  is 


a  s(0.2E- 3)256 


(3.12) 


d_ 

2  fanCT 


(3.13) 


where  d  is  in  unit  of  pixel, 


tana  =  I  d 
2  2  X 


(3.14) 


tana  =  I  tana 
2  25&  2 


(3.15) 


R  =_D  t 
2  fang 
2 


(3.16) 


R  =  1 28D 
1  tan  0.02  * 


(3.17) 


When  the  length  D  of  the  object  is  known,  from  equa¬ 
tion  3.17  we  can  estimate  the  distance  from  lens  to  the 
object  and  is  shown  in  Table  II 


Figure  3.10  Range  Calculation. 


The  distance  between  the  lens  and  the  image  plane 
can  be  adjusted  in  order  to  have  a  clear  image  on  the  film. 
The  camera  has  a  field  of  view  (FOV)  angle  as  shown  in 
Figure  3.9.  The  object  has  to  be  in  the  field  of  view  of 
the  camera.  The  determination  of  the  distance  is  shown  in 
Figure  3.10 

For  simplicity  of  calculation  the  inside  of  the 
camera  is  flipped  to  the  same  side  as  the  object  as  shown  in 
Figure  3.10.  when  the  angle  of  the  FOV  isCK,  1/2  is  the 
half  of  the  full  image  size,  D/2  is  the  half  of  the  dimen¬ 
sion  of  the  object,  X  is  the  distance  from  the  lens  to  the 
image,  and  R  is  the  distance  from  the  lens  to  the  object. 

Assume  that  the  distance  1/2,  distance  Im/2,  and 
angleQ^/2  are  known.  Then,  the  distance  R  can  be  determined 


X  =  JjQQ.tan££ 


I 


2 


2 


(3.9) 


tana'  =  lL 


1  A  V 
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A.  VARIATION  OF  SHIP  SUPERSTRUCTURE  WITH  RANGE 


One  of  the  practical  problem  of  using  the  ship  profile 
is  that  it  is  sensitive  to  range  variations.  Close  ship 
profile  has  more  details  than  the  far  away  ship  profile. 
The  dependency  of  the  geometric  size  of  the  profile  with  the 
range  is  discussed  in  this  section. 

Assume  that  the  object  is  centered  on  the  camera  axis. 
The  field  of  view  of  the  camera,  the  number  of  the  pixel  in 
the  image,  the  size  of  the  image,  and  the  size  of  the  object 
are  known.  Our  problem  is  to  determine  the  distance  between 
the  camera  and  the  object. 

1 .  Background 

The  camera  system  is  similar  to  a  human  eyes  system. 
The  light  reflected  from  the  object  goes  into  the  eyes.  The 
image  of  the  object  falls  on  the  retina,  the  signal  is  sent 
to  the  brain  in  some  electrical  from,  and  the  brain  changes 
it  to  a  from  that  human  can  perceive.  But,  in  the  camera 
system  film  or  senser  are  used  to  pick  up  the  image.  The 
function  of  a  camera  is  shown  in  Figure  3.9 
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It  is  found  by  trials  and  errors  using  the  computer 
program  PARAM  that  in  order  to  retain  maximum  information  of 
the  profile  curve,  the  number  of  knot  positions  required 
should  be  in  the  range  of  25  to  35.  The  number  of  knot 
positions  depends  upon  the  value  of  S  which  must  be  set 
accordingly.  The  appropriate  choice  of  S  to  satisfy  the 
condition  stated  above,  is  important.  For  the  class  of  a 
guided  missile  cruiser(CGN) ,  three  ship  images  at  3 
different  ranges  were  selected.  Then,  run  the  appropriate 
program  for  various  S  factors  to  see  how  the  number  of 
knot(N)  will  vary.  The  results  are  shown  in  Figure  4.6. 
Plots  of  N  vs  S  in  Figure  4.7  through  Figure  4.13  show  that, 
in  most  cases,  the  value  of  N  decrease  quite  rapidly  when 
the  value  of  S  is  in  the  range  of  0  to  100,  and  gradually 
for  S  factor  in  the  range  of  100  to  200,  thereafter,  the 
value  of  N  decreases  very  little.  Obviously,  the  curve  seems 
to  decay  exponentially.  Furthermore,  for  some  classes  of 
ships  changes  are  more  pronounced  than  the  others  which  is 
probably  due  to  the  actual  number  of  knots  present  in  the 
profile.  For  guided  missile  cruiser  ship(CGN)  with  2  lumps, 
the  number  of  knot  positions  required  can  be  33  as  shown  in 
Figure  4.6  We  select  the  factor  S  to  be  about  100. 

The  selection  of  the  factor  S  depends  upon  the 
number  of  the  original  sampling  points.  If  the  factor  S  is 
small,  large  number  of  knots  are  needed.  When  the  factor  S 
is  large,  small  number  of  knots  are  needed.  When  the  number 
of  knot  and  the  Cx  and  Cy  coefficients  are  small,  the 
B-spline  coefficients  Cx  and  Cy  obtained  can  not  be  used  to 
reconstruct  the  curve  close  to  the  original  profile. 


Figure  4.6  Plot  N  vs  S  for  a  CGN 
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3.  The  Output  Cx,Cy  and  Knots  Profiles 

From  the  previous  study  the  factor  S  is  selected  for 
each  of  the  class  of  ship  as  follows 

1.  Destroyer(DD) ,  S  =  20.0 

2.  Container,  S  =  100.0 

3.  Freighter,  S  =  100.0 

4.  Replenishment  oiler,  S  =  20.0 

5.  Tank  landing  ship(LST),  S  =  25.0 

6.  Frigate(FF),  S  =  100.0 

7.  Guided  missile  cruiser (CGN) ,  S  =  125.0 

8.  Guided  missile  destroyer (DDG) ,  S  =  125.0 

The  plot  of  the  B-spline  coef f icients , Cx  and  Cy, and 
X,Y  at  the  positions  of  the  sampling  points  vs  the  Z  param¬ 
eter  are  shown  in  Figure  4.14  through  Figure  4.21. 
Observation  and  comparisons  of  the  curves  show  that  the 
values  of  Cx  and  Cy  exhibit  changes  similar  to  that  of  X  and 
Y  except  that  the  variation  of  values  leads  that  of  the  X 
and  Y.  This  is  due  to  the  fact  that  Cx  and  Cy  have  to  act 
as  the  controlling  factor  for  the  reconstructed  B-spline 
curve  to  get  the  result  close  to  the  original  curve. 

Examination  of  the  plots  of  the  results  of  X  and  Y 
show  that  when  X  is  increasing  monotonically,  Y  is  almost 
constant;  but  when  X  is  almost  constant,  Y  is  increases  or 
decreases.  This  behavior  relating  to  profile  reconstruction 
may  be  explained  as  follows.  For  a  ship  profile  when  X  is 
increasing  and  Y  not  increase  too  much,  this  may  be  inter¬ 
preted  as  an  almost  leveled  profile.  When  X  is  almost 
constant,  and  Y  may  be  increasing  or  decreasing,  it  may  be 
interpreted  as  the  beginning  or  the  ending  of  the  lump. 
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There  are  two  distinct  procedures  used  in  dealing 
with  two  different  kinds  of  lumps  big  or  small.  Thus,  it  is 
necessary  to  distinguish  in  the  first  place  whether  the  lump 
is  big  or  small.  For  the  big  lump,  the  differences  of  Cy  is 
positive  for  all  initial  four  knots.  It  continues  to  the 
peak  and  decreases  toward  the  ending  of  the  lump.  For  the 
small  lump,  the  difference  of  Cy  is  positive  for  the  first 
two  knots  and  stay  constant  or  positive  for  the  third  knot, 
but  the  difference  of  Cy  for  the  forth  knot  is  negative, 
thereafter,  the  program  proceeds  to  establish  the  following 
values  for  each  lump  detected: 

First,  the  knot  positions  (Z  value)  at 

1.  the  beginning  of  the  lump 

2 .  the  ending  of  the  lump 
Second,  the  Knot  number  for 

1 .  the  beginning  of  the  lump 

2.  the  ending  of  the  lump 
Third, Number  of  lumps  detected. 


tions  as  Cy  is  exhibiting  monotonic  increasing  trend.  When 
the  difference  in  Cx  is  decreasing  or  zero,  the  value  of  Cy 
will  have  a  pronounced  change  where  the  beginning  or  the 
ending  points  of  a  lump  can  be  detected.  The  value  of  knot 
position(T)  at  those  points  related  to  the  sizes  of  the 
Cx,Cy  can  be  determined.  Finally,  with  the  values  of  Cx  and 
Cy  at  those  points  known,  the  area  of  the  lumps  can  also  be 
determined. 

Hence,  from  the  ship's  characteristics  and  information 
derived  from  the  above  procedure,  classification  of  ships 
can  be  made  by  considering 

First,  the  number  of  lumps  detected,  1,  2,  or  3 . 

1.  The  number  of  lump=l:  Frigate,  Tank  landing  ship 

2.  The  number  of  lump=2 :  Destroyer,  Guided  missile 
cruiser.  Guide  missile  Destroyer,  and  Replenish¬ 
ment  oiler (AOR) 

3.  The  number  of  lump=3 :  Freighter  and  Container 
Second,  the  position  of  a  lump  relative  the  midships  is 

measured.  This  quantity  is  scaled  by  the  total  length  of 
the  ship.  This  scaled  quantity  will  be  invariant  with 
respect  to  the  different  ship  sizes  at  different  ranges 

Third,  the  area  of  the  lump  is  normalized  to  the  ship 
length  squared. 

1 .  Lump  Detection 

As  shown  in  the  plot  of  X,  Y,  Cx,  and  Cy  vs  Z  param¬ 
eter  in  Figure  4.15,  when  the  difference  (&Cx)  between 
successive  value  of  Cy  varies  from  increasing  to  decreasing, 
and  then,  to  increasing  again,  Cy  exhibits  noticable  varia¬ 
tion.  From  observation  of  Cx  values,  it  is  seen  that  the 
difference  (ACx)  always  has  variation  in  the  same  sequence 
as  stated  above.  Therefore,  only  Cy  values  are  taken  into 
consideration  in  the  program  that  detects  lump. 


7.  Replenishment  oiler  (AOR)  -  There  are  2  little  lumps 
with  the  first  one  starting  at  4/5  of  the  ship  length 
from  the  midships  to  the  left  side  exhibiting  rapid 
decrease  in  height  to  the  level  between  the  bow  and 
the  stern.  The  second  lump  starts  at  1/4  of  the  ship 
length  from  the  midships  to  the  right  side  with  slow 
decline  to  the  point  near  the  stern. 

8.  Freighter-with  3  noticeable  lumps,  the  first  two 
represent  lifting  crane  with  the  first  high  lump 
approximately  1/6  of  the  ship  length  from  the 
midships  to  the  left  side  but  narrow;  meanwhile,  the 
second  lump  is  located  at  approximately  the 
midships  with  the  same  size  as  the  first  one.  The 
beginning  of  the  third  is  at  approximately  1/3  of 
the  ship  length  from  the  midships  to  the  right 
side  with  large  size.  It  is  higher  than  the  first 
two  and  with  gradual  decline  to  the  stern. 

From  the  lump  characteristic  of  different  type  of  ships, 
we  could  distinguish  the  type  of  a  ship  based  on  the  rotated 
superstructure . 

E.  SHIP  CLASSIFICATION 

An  important  aspect  in  categorizing  types  of  ships  is  to 
recognize  the  lumps  above  the  main  deck.  Therefore,  it  is 
useful  to  plot  the  positions  of  X  and  Y  vs  the  Z  parameter 
where  the  Z  parameter  can  be  determined  from  eq(4.5)  where 
Z(I)  is  an  array  (1..M)  as  shown  in  Figure  4.27  where  X  is 
ploted  in  solid  line  and  Y  is  ploted  in  dash  line.  Then, 
plot  the  B-spline  coefficients  Cx  and  Cy  vs  T(N);  the  knots 
position  in  Z  where  N  is  the  total  number  of  knots  as  shown 
in  Figure  4.28.  The  Cx  is  ploted. in  solid  line.  The  Cy  is 
ploted  in  dash  line.  The  Values  of  Cy  will  be  close  to  the 
curve  of  Y  while  the  values  of  Cx  for  the  same  knot  posi- 
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Figure  4-24  Guided  Missile  Destroyer (DDG) . 


Figure  4.25  Destroyer. 


the  first  one.  The  second  lump  ends  at  approximately 
1/4  of  the  ship  length  from  the  midships  to  the  right 
side.  Furthermore,  the  distance  between  the  peak  of 
both  lumps  is  less  than  that  of  the  Guided  missile 
destroyer  shown  .in  Figure  4.26. 
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Figure  4.23  Tank  Landing  Ship(LST). 


Guided  missile  destroyer  -  the  beginning  of  the  lump 
starts  at  approximately  1/4  of  the  ship  length  from 
the  midships  to  the  left  side  with  its  highest  point 
at  the  mast.  After  this  the  slope  will  decline  in  a 
very  rapid  fashion  following  a  noticable  deck 
distance,  then  the  beginning  of  the  second  lump 
occurs  due  to  the  redome  presence.  Therefore,  the 
lump  will  be  narrow  with  great  height,  and  will 
terminate  at  approximatly  1/6  of  the  ship  length  from 
the  midships  to  the  right  side.  Furthermore,  there 
is  a  little  lump  near  the  stern,  this 
distinguish  destroyer  of  the  same  size  as  shown  in 
Figure  4.24. 

Destroyer- lump  characteristic  will  be  similar  to  that 
of  guided  missile  destroyer  except  for  the  small  lump 
as  in  Figure  4.25. 

Guided  missile  cruiser  -  The  beginning  of  the  first 
lump  is  at  1/12  of  the  ship  length  from  the  midships 
to  the  left  side  with  highest  point  at  the  mast.  The 
lump  size  is  large  both  in  length  and  height  and  its 
height  decreases  to  the  point  which  is  a  little  above 
the  level  between  the  bow  and  the  stern.  Therefore, 
the  second  lump  begins  with  almost  the  same  size  as 


D.  SHIP  DESCRIPTION 


The  shape  of  the  ships  depend  upon  the  shapes  and  posi¬ 
tions  of  the  lumps.  Characteristics  of  the  lumps  for 
different  type  of  ships  are  as  follows: 

1.  Frigate  -  The  beginning  of  the  lump  is  at  1/6  of  the 
ship  length  to  the  left  side  of  the  midships  and 
the  peak  of  the  lump  is  at  the  mast  which  is  located 
at  the  midships.  The  termination  of  the  lump 

is  approximate  1/3  of  the  ship  length  from  the 

midships  to  right  side.  In  addition,  the  average 
height  of  the  lump  is  a  little  higher  than  the 

level  between  the  bow  and  the  stern,  and  its  size 
is  1/2  of  the  ship  length  as  shown  in  Figure  4.22. 


Figure  4.22  Frigate. 

2.  Container  -  The  beginning  of  the  lump  is  at  1/3  of 
the  ship  length  to  the  right  side  of  the  midships 
while  the  lump  is  high  and  terminate  at  the  stern. 
The  lump  appears  to  be  in  a  rectangular  shape  with 
small  crane. 

3.  Tank  landing  ship(LST)  -  The  beginning  of  the  lump  is 
at  1/4  of  the  ship  length  from  the  midships  to  the 
left  side;  the  height  of  the  lump  is  higher  than  the 
level  between  the  bow  and  the  stern  by  a  small 
margin,  while  its  highest  point  is  located  at 
approximate  1/6  of  the  ship  length  from  the  midships 
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Plot  X,Y,Cx,  and  Cy  vs  Z  for  a  CGN  at  a  Range  of  45  K-ft. 


a.  Procedure  to  Detect  the  Lump 


1.  Check  a  big  or  small  lump  by  testing  the  conditions 
of  three  varying  values  of  Cy  increments  (ACy)  in 
sequence.  If  they  are  positive  it  is  a  big  lump  and 
set  the  flag- lump  to  1,  otherwise  it  is  a  small  lump 
and  set  the  flag- lump  to  zero;  then  go  to  2 

2.  Check  the  present  knot  position  to  see  where  its 
Z  value  equal  to  zero  or  to  maximum  Z  value,  if  it  is 
Z  maximum  then  stop,  it  not  go  to  3 

3.  If  the  process  begins  at  the  first  knot  position,  set 
the  begin  and  the  flag-end  to  zero;  check  the  status 
of  the  flag- lump  for  1  (big  lump)  or  0  (small  lump), 
if  it  is  a  big  lump,  then  go  to  4.  If  it  is  a  small 
lump ,  then  go  to  7 . 

4.  Check  the  status  for  the  begining  or  ending  of  the 
lump.  If  the  flag-begin  is  1,  it  represents  that  the 
beginning  of  a  lump  is  found,  then  go  to  5.  Otherwise 
it's  not  found,  then  go  to  6. 

5.  Find  the  ending  of  the  big  lump  by  testing  the 

conditions  of  3  values  of  Cy  increments  in  sequence. 

The  first  2  should  be  negative  and  the  third  should 
be  constant  or  positive.  If  the  condition  are 
satisfied,  store  the  Z  value  of  the  position  of  the 
third  knot,  and  the  flag-begin  to  zero;  then  go  to 
10. 

6.  Find  the  beginning  of  the  big  lump  by  testing  the 

conditions  of  3  different  values  of  Cy  increment 

(ACy)  in  sequence.  The  first  Cy  should  be  negative 
or  constant,  the  second  should  be  positive,  and  the 
third  should  be  positive.  If  the  conditions  are 
satisfied,  store  the  Z  value  of  the  position  of  the 
second  knot,  and  set  flag-begin  to  1;  then  go  to  10. 

7.  Check  the  status  for  the  beginning  or  ending  of  the 

lump.  If  the  flag-begin  is  1,  it  represents  the 
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beginning  of  a  lump  is  found,  then  go  to  8.  Otherwise 
it's  not  found,  then  go  to  9. 

8.  Find  the  ending  of  the  small  lump  by  testing  the 

conditions  of  3  values  of  Cy  increment  Cy)  in 
sequence.  The  first  one  should  be  negative  or 

constant.  If  the  conditions  are  satisfied,  store  the 
Z  value  of  the  position  of  the  second  knot,  and  set 
the  flag-begin  to  zero;  then  go  to  10. 

9.  Find  the  beginning  of  the  small  lump  by  testing  the 
conditions  of  3  values  of  Cy  increment  (A  Cy)  in 
sequence.  The  first  should  be  constant  or  negative, 
the  second  should  be  positive,  the  third  should  be 
constant.  If  the  conditions  are  satisfied,  store  the 
position  of  the  second  knot,  and  set  the  flag-begin 
to  1;  then  go  to  10. 

10.  Move  to  the  next  knot,  then  go  to  2. 

This  procedure  is  shown  in  Figure  4.29  and  the  detail  of 
each  procedure  is  shown  in  Appendix  D. 


Figure 


2.  To  Determine  the  Area  Under  a  Lump 

The  area  under  the  lump  can  be  determined  by 


A  AREA  =  ACxACy  +  CyACy 


(4.10) 


Suppose  there  is  a  lump  as  shown  in  Figure  4.30.. 
The  area  increment  can  be  calculated  by  using  eq(4.10).  Then 
the  area  under  BC( which  is  negative)  is  added  to  the  area 
under  AB. 


A  Q 


Figure  4.30  First  Procedure  to  Determine  the  Area. 

Next,  area  under  CD  is  calculated,  which  is  positive 
and  add  this  to  the  last  resulting  area.  Obviously,  the  sum 
would  represent  the  total  area  of  the  lump  as  shown  in 
Figure  4.31. 
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Figure  4.31  Step  by  Step  Procedure  to  Determine  the  Area. 

F.  CONSTRUCTION  OF  THE  DECISION  TREE 

The  characteristics  of  ships  used  for  classification 
are,  the  number  of  lumps,  the  area  of  lumps,  the  knot  posi¬ 
tions^  value)  and  where  the  beginning  of  the  lump  and  lump 
maxima  relative  to  the  midships  are  located. 

In  view  of  the  above  characteristic,  it  is  necessary  to 
find  the  relationships  among  them  to  make  classification 
possible.  Therefore,  for  each  of  the  eight  different  class 
of  ships,  the  decision  tree  is  constructed  which  is  based 
upon  relationships  observed  in  the  plots  of  the  knot  posi¬ 
tion  (Z  values)  for  the  beginning  of  the  lumps  normalized  by 
the  total  ship  length  (Z  value)  VS.  the  number  of  lumps;  the 
area  of  the  lumps  normalized  by  the  total  ship  length 


(2  value)  squared  vs  the  number  of  lumps;  and  the  Z  value  of 
the  lump  maxima,  as  shown  in  Figure  4.32  through 
Figure  4.39. 

Thus,  according  to  the  number  of  lumps  presented  in  the 
profile  used  as  the  first  criteria,  ships  can  be  divided 
into  3  distinct  groups.  Then,  classification,  can  be  made 
by  comparing  further  characteristic  as  shown  in  Table  III, 
and  the  complete  decision  tree  constructed  from  Table  III  is 
shown  in  Figure  4.40. 
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0.00  5 


•  i 
2 


I  2  3  m 


AREA  OF  LUMP 


/  =CONTAINER  AT  RANGE  28C 
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END 


G.  SUMMARY 


The  decision  tree  which  is  constructed  from  Table  III 
does  not  provide  100-percent  correct  classification  for  all 
types  of  ship  images.  Furthermore,  the  presence  of  noise  in 
images  may  cause  complicatations  in  classifying  the  ships. 
For  example,  with  excessive  noise  in  the  original  image, 
Sobel  operator  for  edge  enhancement  in  the  preprocessing 
process,  still  yield  result  with  residual  noise  present  in 
Figure  4.41.  These  residual  noise  is  undesirable  since  it 
causes  failure  to  extract  profiles  which  retain  necessary 
informations  from  the  original  images.  The  subsequent  clas¬ 
sification  of  profile  by  the  Fourier  Coefficient  method  and 
the  B-spline  Coefficient  method  becomes  difficult. 
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Figure  4.41  Noisy  Image. 

As  discussed  before  in  order  to  reduce  the  noise,  an 
appropriate  threshold  gray  value  for  the  preprocessed  image 
is  set.  This  results  in  a  silhouette  image.  However,  if  the 
value  is  too  high  the  profile  becomes  broken  which  prevent 
successful  operation  in  the  closing  process  discussed  in 
chapter  2,  as  shown  in  Figure  4.42.  On  the  other  hand 
decreasing  the  threshold  value  results  in  erroneous  profile, 
as  shown  in  Figure  4.43.  In  some  cases,  when  the  closing 
process  is  used,  certain  vital  information  is  lost.  This 
effect  can  be  seen  in  comparing  the  image  in  Figure  ‘.44, 


and  Figure  4.45.  Therefore,  loss  of  informations  due  to  the 
attempt  of  eliminating  noise  and  the  inability  in  the 
preprocessing  process  to  cope  with  the  residual  noise,  yield 
erroneous  profile.  Consequently,  failure  occurs  in  the 
succeeding  classification  steps. 


Figure  4.45  After  Closing 


V.  CONCLUSION 


The  shape  of  the  superstructure  of  a  ship  is  the  most 
important  feature  used  in  the  classification  algorithms.  To 
extract  the  edge  profile,  a  Sobel  operator  is  employed.  The 
limitation  of  a  Sobel  operator  is  that  some  small  details  of 
the  edge  is  lost.  For  example,  the  image  of  a  DD  at  a  range 
of  79000  feet  after  applying  the  Sobel  operator  shows  the 
superstructure  and  the  radar.  But  the  edge  of  a  small  mast 
disappeared  as  shown  in  Figure  5.1.  Furthermore,  if  the 
threshold  value  is  set  too  big,  the  edge  image  of  the  super¬ 
structure  profile  becomes  broken.  The  top  superstructure 
profile  is  obtained  by  setting  the  gray  value  under  the 
slope  between  the  bow  and  the  stern  to  zero.  We  need  to 
apply  a  contour  tracking  process  to  refine  the  superstruc¬ 
ture  profile.  For  some  images,  the  connection  of  the  broken 
profile  pieces  may  be  achieved  in  a  Closing  operation  as 
discussed  in  Chapter  2.  The  disadvantage  of  the  Closing 
operation  is  that  some  small  details  of  the  profile  may 
disappear.  The  superstructure  profile  of  a  freighter  at  a 
range  of  53000  feet  is  shown  in  Figure  5.2.  After  the 
Closing  process  the  detail  of  the  cranes  disappeared  as 
shown  in  Figure  5.3. 

The  ship  classification  can  be  achieved  by  either  the 
Fourier  Coefficient  method  or  the  B-spline  Coefficient 
method.  In  using  these  coefficients  to  classify  ships,  it 
is  found  that  only  the  initial  coefficients  that  lie  between 
the  0th  to  the  20th,  are  relevant,  while  the  rest  are  not. 
Inspection  of  the  comparison  curves  of  the  same  class  of 
ships  shows  that,  similarities  in  patterns  exists  up  to  the 
20-th  point.  Beyond  that,  diversities  in  shape  are  so  great 
that  inclusion  of  those  additional  points  for  classification 


Figure  5.5  Logarithmic  Magnitude  of  a  CGN  at  a  Range  of  55000  feet. 


In  the  B- spline  coefficient  method  the  technique  of 
uneven  knot  selection  has  been  employed.  With  the  original 
sampling  points  on  the  ship  profile  as  input,  a  smaller  set 
of  approximation  samples  are  collected.  These  can  be  used  to 
reconstruct  the  ship  profile  with  sufficient  information 
retained  from  the  original  image.  The  execution  time  as 
approached  to  that  of  handling  the  original  data  directly 
has  been  greatly  reduced.  The  reduction  of  the  number  of 
sample  points  by  almost  a  factor  of  10  is  common.  For  a  CGN 
ship  a  set  of  original  sampling  points  of  290,  has  been 
reduced  to  36. 

Comparing  the  two  methods,  Fourier  Coefficient  and 
B-spline  Coefficient,  the  former  method,  in  some  cases  is 
not  effective  to  establish  satisfactory  classification  of 
ships.  This  is  due  to  difficulties  in  matching  similarities 
of  the  shape  of  the  coefficient  curves.  The  latter, 
however,  surpasses  the  former  in  that  it  is  able  to  classify 
more  ships  accurately  using  computer  programs.  It  is 
possible  to  improve  the  reliability  of  those  two  methods  by 
reducing  noise  in  the  data  collection  process  and  the 
preprocessing  process. 


6-Dec-1984  17:18:40  VAX- 
Source  Listing  6“0ec-1984  17:18:31  _ORA 

program  cutUnput,outoutiinfil«,outfil<)! 

ty  pa 

Oyt*  *  0..255: 

liaagaroail  =  oackad  array  CO. .2573  of  byte; 

iiaag.»ro*2  *  packed  array  CO. .2553  of  integer; 

inageroa.3  -  oacktd  array  CO. .2553  of  byte; 

roal  •  packed  array  CO. .2553  of  byte: 
ro«2  *  oackad  array  C0..2S53  of  integer! 

var 

sooel  :  array  CO. .633  of  iaagerow2 : 
i.J  :  byta; 

f  :  array  CO. .653  of  inagaroal i 
outfile  :  file  of  inagaroal! 
da ■ dy < range. or ignt ■ n ,n  linteger; 
infila  :  file  of  roal! 
slope  :  real: 

iaaga  ’array  CO. .633  of  roal! 
a . al »a2 » yl » y2  :  integer; 

NUMl,NUM2»NUM3,NUM4,iea*,min,«aal  :  integer; 
cal  Sarray  CO. .633  of  roal; 
thas  :  integer; 

name  :  paCkeo  array  CI..203  of  char: 

BEGIN 

RRITELNC "INPUT  SHIP  FILENAME  OUT  CUT. OAT'); 

REAOLNCNAME) 5 
URITcLNC 'THRESHOLD') ! 

'AOCTMES): 

'ELNC'NUH  TC  CUT  FROM  L  cFT ' ) ! 

Rl  iCNUMl); 

HR  I i EL  N( ' N JM2  TO  CUT  FROM  RIGHT*); 

REA0CNUM2): 

HR1TJLNC'NUM3  TO  CUT  FROM  TOP'); 

REA0CNUM3) : 

MRITELNC *NUM4  TO  CUT  FROM  BUTTON'); 

REA0CNUM4) ; 

open  ( in f il e .NAME , hi story  :«old» 
access. aa  t  hod  ^sequential. 

record.length  :* 256. record _ t yp e  :*fi«ed); 
open  ( out  file . 'CUT • dat' .History  : *nei«»record_lenoth  :*256. 

record.tyoe  :*fiaed): 
resetlinfile)  ; 
reerite  (outfile); 
i :  *0 ; 

ahile  not  eof  (infile)  do 
oe  gi  n 

read  ( in f 1  la , leageC i 3) ; 
for  j:*0  to  255  do 

fCi*l.j*13  :«  iaageCi.J3: 
i : «i ♦  1 : 
end: 

(dcoapute  sooele) 
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TEKM2  =  TERM2«CY(L0)sCC:t, J) 

170  CONTINUE 

TERM  =  wC  IT  )*(  (TERH1-X  (I  T)  )o*2*(  7E  RM2-T C  I  T  )  )*S2  ) 

FPART  =  FPART+TERM 
IFCNEW.E3.0)  GO  TO  180 
STORE  *  TERM*0.5 
FPInTCI)  *  FPART-STORE 
I  =  1*1 
FPART  *  STORE 
NEW  =  0 

180  CONTINUE 

F  P I  NT ( NR  I N  T )  •  FPART 
00  190  L  =  1 tNPLUS 
C  A  00  A  NEW  KNOT. 

CALL  NKN3TC2,HfT,N,FPINT.NR0ATA,NRINT) 

C  TEST  WHETHER  WE  CANNOT  FURTHER  INCREASE  THE  NUM9ER  OF  KNOTS. 
1FCN.E3.NHAX  .OR.  N.EQ.NEST)  GO  TO  200 
190  CONTINUE 

C  RESTART  THE  CONFUTATIONS  WITH  THE  NEW  SET  OF  KNOTS. 

200  CONTINUE 

C  TEST  WHETHER  the  APPR OX  I H A T I  ON  X  »  SX  <  Z  )  .  Y  *  5  Y  C  Z  )  WITH  SXCZ3  AND  STCZ) 

C  THE  LEAST-S3UARES  KTH  OcGREE  POLYNOMIALS.  IS  A  SOLUTION  OF  OUR  PROBLEM 
250  IFCIER.E3.-2)  GO  TO  440 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C  PART  2:  DETERMINATION  OF  THE  SMOOTHING  SPLINES  SXPCZ)  AND  SYPCZ3.  C 
C  ***»*$***««  3$  :esi#si*3*sea***se*#*$3*s#:»*3iPse  ****#$::#  »*«**»«»*:  C 
C  WE  HA  VE  OETERMlNEO  THE  NUMBER  OF  KNOTS  AND  THEIR  POSITION.  C 
C  WE  NOW  COMPUTE  The  8-SPLlNE  COEFFICIENTS  OF  THE  SMOOTHING  SPLINES  C 


C  SXPCZ)  AND  SYP(Z).  THE  OBSERVATION  MATRIX  A  IS  EXTENDED  6T  THE  ROWS  C 
C  OF  MATRIX  B  EXPRESSING  THAT  THE  KTH  DERIVATIVE  DISCONTINUITIES  OF  C 
C  SxP(Z)  AND  STPCZJ  AT  THE  INTERIOR  KNOTS  T(K*2),...UN-x-1)  MUST  BE  C 
C  ZERO.  THE  CORRESPONDING  WEIGHTS  OF  THESE  ADDITIONAL  ROWS  ARE  SET  C 
C  TO  1/SQRTCP}.  ITERATIVELY  WE  THEN  HAVE  TO  OETERMlNE  THE  VALUE  OF  P  C 
C  SUCH  THAT  FtP)*SUMCWl*((XI-SXP<ZI))«*2«(ri-SYP(ZI))»F2)  BE  =  S.  WE  C 
C  ALREADY  know  THAT  THE  LEAST- SQUARE S  POLYNOMIALS  CORRESPOND  TO  P*0,  C 
c  AND  that  The  LEAST-SQUARES  SPLINES  CORRESPOND  TO  P=INFINITT.  THE  C 
C  ITERATION  PROCESS  WHICH  IS  PROPOSED  HERE,  MAKES  USE  OF  RATIONAL  C 

C  INTERPOLATION.  SInC!  PCP1  IS  A  CONVEX  ANO  STRICTLY  DECREASING  C 

C  FUNCTION  OF  P,  IT  CAN  BE  APPROXIMATED  BY  A  RATIONAL  FUNCTION  C 

c  rcpj  *  cu*p»v j/<p»w).  three  values  of  p<pi,p2,p3>  with  correspond-  c 
C  ING  VALUES  of  F(P)  <Fl*F(Pn-S,F2=F(P2)-S,F3«F(P3j_s)  ARE  USED  C 

C  TO  CALCULATE  THE  NEW  VALUE  OF  P  SUCH  THAT  R(P)*S.  CONVERGENCE  IS  C 

c  guaranteed  by  taking  fi>o  and  F3<o.  c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  EVALUATE  THE  DISCONTINUITY  JUMP  OF  THE  KTH  DERIVATIVE  Oc  THE 
C  B-SPUNES  AT  THE  KNOTS  T  C  L  3  ,  L  *K  *2  .  .  .  .  N-K- 1  AND  STORE  IN  8. 

CALL  C  ISCOC T ,N.K2 ,B) 

C  INITIAL  VALUE  FOR  P. 

PI  ■  0. 

FI  *  FPO-S 
P3  *  -1. 

F  3  «  F  Pm  S 
*  *  -F1/F3 
ICHECK  *  0 
N8  «  N-NMIM 

C  ITERATION  PROCESS  TO  FINO  THE  ROOT  OF  F<P)  *  S. 

00  350  ITER*1,MAXIT 

C  the  ROWS  of  MATRIX  ?  WITH  WEIGHT  1/SQRT(P>  ARE  ROTATED  INTO 
C  The  TRIANGULARISEO  OBSERVATION  matrix  A  WHICH  is  stored  in  g. 

P 1 NV  *  1.0/P 

00  260  1*1, NK1 


39  30 
3990 
4000 
4010 
4020 
4030 
4040 
4050 
4060 
4070 
4080 
4090 
4100 
4110 
4120 
4130 
4  14  0 
4  150 
4160 
4  170 
4  1  80 
4  1  90 
4200 
4210 
4220 
4230 
4240 
4250 
4260 

42  70 
4  2  B  0 
4290 
4300 
4310 
4320 
4330 
4340 

43  50 
4360 
4370 
4380 
4390 

44  30 
4410 
4420 
4430 
444  0 
4450 
4460 

44  70 
4480 
4490 
4500 
4510 
4  520 
4530 
4540 

45  50 
4560 
4570 
4  c  »r 
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IFCPIV.EQ.0.3  GO  TO  110  3370 

C  CALCULATE  THE  PARAMETERS  OP  THE  GIVENS  TRANSFORMATION,  3380 

CALL  C0SSINtPIV.Wl,0lAGCJ3.C0S.SIN3  3390 

C  TRANSFORMATIONS  to  RIGHT  HAND  SIDES.  3*00 

CALL  ROTATECPIV,COS,SIN,XI,RXCJ33  3*10 

CALL  R3TATE(PIV,C0S,SIH,ri,RT<J))  3*20 

IFCI.EQ.K13  GO  TO  120  3430 

12  *  0  3**0 

13  *  1*1  3450 

00  100  II  *  13. XI  3450 

12  *  12*1  3*70 

C  TRANSFORMATIONS  TO  LEFT  HAND  SIDE.  3490 

CALL  R0TATECPIV.C0S.SIN.HCI13.ACJ.I233  3*90 

100  CONTINUE  3500 

110  CONTINUE  3510 

C  AOD  CONTRIBUTION  OF  THIS  ROW  TO  THE  SUM  OF  SQUARES  OF  RESIOUAL  3520 

C  RIGHT  MANO  SIOES.  3530 

120  FP  *  FP.uI*<Xl*P2*TI»*2)  35*0 

130  CONTINUE  3550 

IFCIER.EQ.-23  FPO  *  FP  3550 

C  BACKWARD  SUBSTITUTION  TO  OBTAIN  THE  B-SPLINE  COEFFICIENTS.  3570 

CALL  BACXCA, RX.NK1.K.CX)  3580 

CALL  BACKC A.RT.NKl.K.CT)  3590 

C  TEST  WHETHER  THE  APPROXIMATION  X * SX I NF C Z 3 . T* S T INF ( 2 3  IS  AN  3600 

C  ACCEPTABLE  SOLUTION.  3610 

FPMS  *  FP-S  3630 

IFCABSCFPMSJ.LT. ACC3  GO  TO  **0  3630 

C  IF  F ( P* INF  3  <  S  ACCEPT  THE  CHOICE  OF  KNOTS.  36*0 

IFCFPMS.LT. 0.3  GO  TO  250  3650 

C  IF  N-NMAX.SXINFCI3  ANO  SYINFCZ3  ARE  INTERPOLATING  SPLINES.  3660 

IFCN.E0.NMAX3  GO  TO  *30  3670 

c  Increase  ThE  number  of  xnots.  3680 

C  IF  N*NEST  WE  CANNOT  INCREASE  THE  NUMBER  OF  KNOTS  BECAUSE  OF  3690 

C  THE  STORAGE  CAPACITY  LIMITATION.  3700 

IFCN.EG.NEST3  GO  TO  *20  3710 

C  OETEPMINE  THE  number  OF  XNOTS  NPLUS  WE  ARE  GOING  TO  AOD.  3T20 

1FCIER.EQ.03  GO  TO  1*0  3730 

NPLUS  *  1  37*0 

IER  *  0  3750 

GO  TO  150  3760 

1*0  NPL1  •  NPLJS*2  3770 

IFCFP0L0-FP.GT.ACC3  NPL1  *  FL 0 A TC NPLUS JPFPHS/ C F POL D-F P )  3780 

NPLUS  *  MIN0CNPLUS*2.MAX0CNPL1.NPLUS/2,1)3  3790 

ISO  FPOLO  *  FP  3800 

C  COMPUTE  THE  SUMCWI*CCXI-SXINFCZI33**2»CYI-SYINFCZI)3**233  FOR  3810 

C  EACH  XNQT  INTERVAL  TCJ«K3  <*  21  <*  TCJ*K»13  AND  STORE  IT  IN  3820 

C  FP1NTCJ3.J*1.2....NRINT,  3830 

FPART  *  0.  38*0 

I  .  i  3850 

L  *  K2  3860 

NEW  »  0  3870 

DO  180  I T  *  1 , M  3880 

IFCICIT).LT.TCL3  .OR.  L.GT.NK1)  GO  TO  160  3890 

NEW  *  1  3900 

L  *  L ♦ l  3910 

160  TERM!  *  0.  3920 

TERM2  =  0.  3930 

LO  *  L-K2  39*0 

00  170  J»l,Kl  3950 

LO  *  L  0» 1  3960 

T  E  RMl  s  7ERM1*CX(L03PCC:T.J3  ’'*Tr 
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40  CONTINUE  2750 

CO  TO  50  2770 

C  IF  S>0  OUR  INITIAL  CHOICE  OF  KNOTS  CEPENDS  ON  THE  VALUE  OF  IOPT.  2790 

C  IF  IOPT *0  OR  IOPT* I  ANO  S>*FP0,  WE  START  COMPUTING  THE  L E AS T - SOuA R E S  2790 

C  POLYNOMIALS  Oc  DEGREE  K  WHICH  ARE  SPLINES  WITHOUT  INTERIOR  KNOTS.  2800 

C  IF  I0PT*1  ANO  FPO>S  WE  START  COMPUTING  THE  L E A S T - SOU ARE S  SPLINES  2810 

C  ACCORDING  TO  THE  SET  OF  KNOTS  FOUNO  AT  THE  LAST  CALL  OF  THE  ROUTINE.  2820 

45  IFCIOPT.LE.O)  GO  TO  SO  2830 

IFCFPO.GT.S)  GO  TO  60  2840 

SO  N  *  NMIN  2850 

NROATA(l)  *  M-2  2860 

C  MAIN  LOOP  FOR  THE  DIFFERENT  SETS  OF  KNOTS.  M  IS  A  SAVE  UPPER  BOUND  2870 

C  FOR  THE  NUMBER  OF  TRIALS.  2890 

60  00  200  ITER  *  1 « M  2890 

IFCN.E3.NMIN)  IER  *  -2  2900 

C  FINO  NRINT.  TNE  NUM3ER  OF  KNOT  INTERVALS.  2910 

NRINT  *  N-NMIN«1  2920 

C  FINO  THE  POSITION  OF  THE  ADDITIONAL  KNOTS  WHICH  ARE  NEEDED  FOR  2930 

C  THE  8-SPLINE  REPRESENTATION  OF  SXCZ)  ANO  STC2).  2940 

NKI  *  N-Kl  2950 

I  *  N  2960 

00  70  J*1.KI  2970 

TCJ)  *  ZB  2980 

TCI)  *  ZE  2990 

I  *  1-1  3000 

70  CONTINUE  3010 

C  COMPUTE  THE  B-SPLINE  COEFFICIENTS  OF  THE  LE AST-SOUARES  SPLINES  SXINFCZ)  3020 
C  ANO  STINFCZ).  THE  OBSERVATION  MATRIX  A  IS  BUILT  UP  ROW  3T  ROW  AND  3030 

C  RSOUCEO  TO  UPPER  TRIANGULAR  FORM  9Y  GIVENS  TRANSFORMATIONS  3040 

C  WITHOUT  SQUARE  ROOTS.  AT  THE  SAME  TIME  FP*FCP»INF)  IS  COMPUTED  3050 

FP  *  0.  3060 

C  INITIALIZE  The  OBSERVATION  MATRIX  A.  3070 

03  80  I«1,NK1  3090 

OIAGCI)  *  0.  3090 

RXCI)  *  0.  3100 

RTCI)  ■  0.  3110 

00  80  J*1,K  3120 

ACI.J)  *  0.  3130 

80  CONTINUE  3140 

L  ■  K1  3150 

00  130  17*1, M  3160 

C  FETCH  THE  CURRENT  DATA  POINT  X C I T ) , T C I T )  ,  Z C  I T )  .  3170 

XI  »  XCIT)  3180 

T I  *  TCIT)  3190 

Z1  ■  ZCIT)  3200 

WI  »  WCIT)  3210 

C  SEARCH  FOR  KNOT  INTERVAL  TCL)  <*  ZI  <*  TCL+1).  3220 

IFCZI.GE.TCL*! )  .AND.  L.NE.NK1)  L  *  L«1  3230 

C  EVALUATE  THE  C K ♦ 1 )  NON-ZERO  8-SPLINES  AT  Z I  AND  STORE  THEM  IN  Q.  3240 

CALL  BSPLINCT.N.K.ZI.L.M)  3250 

00  90  1*1, K1  3260 

IFCMCD.LT.0.1E-07)  MCI)  «  0.  3270 

QCIT.I)  «  MCI)  3290 

90  CONTINUE  3290 

C  ROTATE  THE  NEW  ROW  OF  THE  OBSERVATION  MATRIX  INTO  TRIANGLE  BT  3300 

C  GIVENS  TRANSFORMATIONS  WITHOUT  SQUARE  ROOTS.  3310 

J  *  t-Kl  3320 

00  110  1*1, Kl  3330 

IFCWI.EO.O.)  GO  TO  130  3340 

J  *  J*l  3350 

PI  V  *  MCI)  ’ 3o0 
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IPCK.Ce-03  IER  *  10 

IFCH.LT.IU  .OR.  NEST.LT.NH1N)  IER  *  10 
IFCS.LT.O.)  IER  *  10 
IFCIER.NE.O)  CO  TO  *4Q 

C  CHECK  WHETHER  THE  Z-VALUES  ARE  PROVIDED  WITH  BT  THE  USER. 
IFCIPAR.NE.O)  CD  TO  6 

C  FINO  FOR  EACH  OATA  POINT  A  CORRESPONDING  VALUE  OF  THE  PARAMETER  Z 
C  ANO  FIX  THE  BOUNDARIES  ZB  ANO  ZE. 

Ul)  «  0. 

00  4  1*2, H 

ZCI)  «  ZCI-1)*SQRTCCXCI)-XCI-1)>**2*CVCI>-YCI-1>)*P2) 

4  CONTINUE 
ZB  *  ZCI) 

ZE  *  ZCH) 

6  IFCZ8.GT.ZC1)  .OR.  ZE.LT.ZCM)  .OR.  WCD.LE.O.)  IER  *  10 

00  10  1*2, H 

IFCZa-D.CE.ZCI)  .OR.  WCD.LE.O.)  IER  *  10 
10  CONTINUE 

IFCIER.NE.O)  GO  TO  440 

C  CALCULATION  OF  ACC,  THE  ABSOLUTE  TOLERANCE  FDR  THE  ROOT  OF  FCP)=S. 

ACC  *  TOL*S 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C  PART  l:  DETERMINATION  OF  THE  NUMBER  OF  KNOTS  ANO  THEIR  POSITION  C 
C  »M*****<«*******9*******«*>«***«««»*«P****«*P**«*«9*««»«**«t«  C 
C  GIVEN  A  SET  OF  KNOTS  WE  COMPUTE  THE  LE A ST-S3U AR E S  SPLINES  SXINFCZ)  C 
C  ANO  STINFCZ).IF  THE  SUM  FCP*INFX*S  WE  ACCEPT  THE  CHOICE  OF  KNOTS.  C 
C  OTHERWISE  WE  HAVE  TO  INCREASE  ThEIR  NUMBER.  C 
C  THE  INITIAL  CHOICE  OF  KNOTS  OEPENOS  ON  THE  VALUE  DF  S  ANO  IOPT.  C 
C  IF  S*0  WE  HAVE  SPLINE  INTERPOLATION!  IN  THAT  CASE  THE  NUMBER  OF  C 
C  KNOTS  E3UALS  NMAX  »  M*K*1.  C 
C  IF  S  >  0  ANO  C 
C  IOPT *0  WE  FIRST  COMPUTE  THE  LEAST-SQUARES  POLYNOMIALS  OF  C 
C  OEGREE  K!  N  »  NMIN  «  2*K«2  C 
C  I OPT* 1  WE  START  WITH  THE  SET  OF  KNOTS  FOUND  AT  THE  LAST  C 
C  CALL  OF  THE  ROUTINE,  EXCEPT  FOR  THE  CASE  THAT  S  >  FPO!  THEN  C 
C  WE  COMPUTE  DIRECTLY  THE  LEAST-SQUARES  POLYNOMIALS  OF  DEGREE  K.  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C  OETERMINE  NMAX,  THE  NUMBER  OF  KNOTS  FOR  SPLINE  INTERPOLATION. 

NMAX  *  M»K1 
IFCS.GT.O.)  GO  TO  45 

C  IF  S«0.  SXCZ)  AND  SYCZ)  ARE  INTERPOLATING  SPLINES. 

N  *  NMAX 

C  TEST  WHETHER  THE  REQUIRED  STORAGE  SPACE  EXCEEDS  THE  AVAILABLE  ONE. 
IFCN.GT.NEST)  GO  TO  420 

C  FINO  THE  POSITION  OF  THE  INTERIOR  KNOTS  IN  CASE  OF  INTERPOLATION. 

MK1  ■  M-Kl 

IFCHK1.EQ.0)  GO  TO  60 
K 3  *  K/2 
I  ■  K2 
J  •  K3»2 

IFCK3*2.EC.K )  GO  TO  30 
DO  20  L*1,MK1 
TCI)  *  ZCJ) 

I  «  1*1 
J  *  I 
20  CONTINUE 
GO  TO  60 

30  DO  40  L • 1 , MK 1 

TCI)  *  CZCJ)-ZCJ-1))*0.5 
I  *  1*1 
)  «  J«1 


21  SO 
2160 
2170 
2130 
2190 
2200 
2210 
2220 
2230 
2240 

22  SO 
2260 
2270 
2280 
2290 
2300 
2310 
2320 
2330 
2340 
2350 
2360 
2370 
2380 
2390 
2400 
2*10 
2420 
2430 
2440 
2450 
2460 
2470 
2480 
2490 
2500 
2510 
2520 
2530 
2540 
2550 
2560 
2  5  T  0 
2580 
2S90 
2600 
2610 
2620 
2630 
2640 
2650 
2660 
2670 
2680 
2690 
2700 
2710 
2  720 
2730 
2740 

,  T  <  p 
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•  —’Ti  w  »  'Ti  V  V 
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IER»-2 : NORMAL  RETURN,  SXCZ)  ANO  ST(2)  ARE  POLYNOMIALS  OP  OEGREE  K 
IE  R  >0  lABNORMAL  TERMINATION. 

IER*i:THE  REQUIRED  STORAGE  SPACE  EXCEEDS  THE  AVAILABLE 
STORAGE  SPACE, SPECIFIED  THE  PARAMETER  NEST. 

PROBABLY  CAUSES! NEST  OR  S  TOO  SMALL. 

IER*2:A  ThEORETICALLT  IMPOSSIBLE  RESULT  WAS  FOUND  DURING 
THE  ITERATION  PROCESS. 

PROBABLY  CAUSES :  TOL  TOO  SMALL 
IER*3! The  MAXIMAL  HUMBER  OF  ITERATIONS  MAXIT  MAS  BEEN 
REACHED. 

PROBABLY  CAUSES! MAXIT  OR  TOL  TOO  SMALL. 

IER* 1 0  ?  SOME  OF  THE  INPUT  OATA  ARE  INVALIDCSEE  RESTRICTIONS). 

RESTRICTIONS: 

1)  M  >  K  >  0 

2)  IB  <*  ICR)  <  I< R* 1 )  <«  IE,  R*1,2,...H-1.  CIPAR  *  1) 

3)  NCR)  >  0,  R« 1 , 2 , .. • M. 

A)  S  >■  0. 

5)  NEST  >*  2*K*  2. 

other  subroutines  required: 

BSPLlN.COSSIN, ROTATE, BACK. NXNOT, DISCO  AND  RATION. 

DIMENSION  XCM),Y<N),U<M),2(M),TC20Q),CX(2C0).CYC200), 

<  FPINT(200).RX<200),RVC200).QIAGC200>,OPAIM£C200), 

<  GC20C,6).BC200,7),QC*00,6),M{7) .NR0ATAC200) , A (200, 5) 

COMMON/OPT 1/NRO AT ACNE ST ), FPO ,F POLO. NPLUS 

nroata:  INTEGER  ARRAY.LENGTH  nest, WHICH  GIVES  THE  NJMBER  of 
OATA  POINTS  INSIOc  EACH  KNOT  INTERVAL. 

FPO  !  REAL  VALUE,  WHICH  CONTAINS  THE  SUMC W I  * < « I - S X C Z I > )*»2 ) ♦ 

SUM(WI*CTI-SYCII))*P2)  WITH  SXCI)  ANO  SYCI)  LE AST- SOU AR E S 
POLYNOMIALS  OF  OEGREE  X. 

F POLO  5  REAL  VALUE, WHICH  CONTAINS  THE  SUM< WIP C X I- S X C I I ) >**Z >♦ 

SUM(wI»CYI-SYC2I))**2)  with  SXC2)  ANO  SYCI)  LEAST-SauARES 
SPLINE  FUNCTIONS  CORRESPONDING  TO  THE  LAST  FOUND  SET  DF 
KNOTS  BUT  ONE. 

NPLUS  :  INTEGER  VALUE. GIVING  THE  NUMBER  OF  KNOTS  OF  THE  LAST 
SET  MINUS  THE  NUMBER  OF  THE  LAST  SET  BUT  ONE. 
COMMON/OPTI/NROATA.FPO.FPOLO.NPLUS 
OATA  INITIALIZATION  STATEMENT  TO  SPECIFY 

TOL  !  THE  REQUESTED  RELATIVE  ACCURACY  FOR  THE  ROOT  Oc  FCP)  »  S. 
MAXIT:  THE  MAXIMAL  NUMBER  OF  ITERATIONS  ALLOWED. 

NEST  :  AN  OVER-ESTIMATE  OF  THE  NUMBER  OF  KNOTS  N.  THIS  PARAMETER 
must  BE  SET  BY  THE  USER  TO  INOICATE  THE  STORAGE  SPACE 
AVAILABLE  TO  THE  SUBROUTINE,  the  DIMENSION  SPECIFICATIONS 
OF  THE  ARRAYS  T , C X , C Y , NRO A T A , FPINT , RX , R Y , DI AG , OPRI ME C N) , 
A(N,K),G(N,K«1),6(N,K*2>,C(M,K«1)  AND  H ( K ♦ 2 )  DEPEND 
ON  N,M  ANO  K.  SINCE  N  IS  UNKNOWN  AT  THE  TIME  THE 
USER  SETS  UP  THE  DIMENSION  INFORMATION  AN  OVER-ESTIMATE 
OF  THESE  ARRAYS  WILL  GENERALLY  BE  MADE.  THE  FOLLOWING 
REMARKS  ARE  INTENDED  TO  HELP  THE  USER 

1)  2*K«  2  <«  N  <*  M4K ♦ I 

2)  THE  SMALLER  THE  VALUE  OF  S,  THE  GREATER  N  WILL  BE. 

3)  NORMALLY  N  *  M/2  IS  AN  OVER-ESTIMATE. 

OATA  T0L/0.001/.NAXIT/20/.NEST/200/ 

BEFORE  STARTING  COMPUTATIONS  A  OATA  CHECK  IS  MADE.  IF  THE  INPUT 
OATA  ARE  INVALIO  CONTROLE  IS  IMMEDIATELY  REPASSED  TO  THE  DRIVER 
PROGRAM  <IER*10). 

IER  ■  0 
K  1  •  K  ♦  I 
K  2  «  K  1  ♦  1 
NMIN  «  2*K 1 


1540 
1550 
1560 
1570 
1  5  80 
1590 
1600 
1610 
1620 
1630 
1640 
1  6S0 
1660 
1670 
1 6B0 
1690 
1700 
1710 
1720 
1730 
1740 
1750 


1790 
1600 
1810 
1820 
1 B  30 
1B40 
1650 
I860 
1870 
1880 
1890 
1900 
1910 
1920 
1930 
1  940 
t  950 
1  960 
1970 
1980 
1990 
2000 
2010 
2020 
2030 
2040 
2050 
2060 

2080 
2090 
2100 
2110 
2120 
2130 
?  1  AC 


-H 


■  •) 
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SUBROUTINE  PARAHCX,T,Z,W,M,ZB,ZE,K,S.N,T,CX,CY.FP,I0PT,IPAR,IER)  0930 

GIVEN  THE  SET  OF  OATA  POINTS  (KIJ.KD)  W ITH  CORRESPONDING  2-  09*0 

VALUES  ZCI) ,1*1. 2. ...M  ANO  GIVEN  ALSO  THE  SET  OF  POSITIVE  0950 

NUMBERS  WCI),I*1.2....H,  SUBROUTINE  PARAM  FINOS  A  SMOOTH  APPROXIMA1-  0960 

ING  CURVE  WITH  PARAMETER  REPRESENTATION  X  *  SX(Z>.  T  =  SYC2).  0970 

SXCZ)  ANO  STCZ)  ARE  TWO  SPLINE  FUNCTIONS  OF  DEGREE  X  WITH  THE  NUMBER  0980 

ANO  THE  POSITION  OF  THE  KNOTS  TCJ), J-1.2....N  AUTOMATICALLY  0990 

chosen  by  the  routine,  the  smoothness  of  sxcz>  and  sycz>  is  1000 

ACHIEVED  BY  MINIHALIZINS  THE  SUMC OX C R )*#2 *0 Y C R ) PP2 )  WHERE  DXCR)  1010 

ANO  DVCR)  STAND  FOR  THE  OISCONT I NU I T Y JUMP  OF  THE  XTH  DERIVATIVE  1020 

OF  SXCZ)  ANO  STCZ)  AT  THE  KNOT  TCR).  R«K *2  , . . . N-K-l .  1030 

THE  AMOUNT  OF  SMOOTHNESS  IS  DETERMINED  BY  THE  CONDITION  THAT  FCP)  «  10*0 

SUMCWCI)*CCXCI)-SXCZCI)))R*2  ♦CYCI)-STCZCI)))**2))  BE  <*  S,  WITH  1050 

S  A  GIVEN  NON-NEGATIVE  CONSTANT.  1060 

THE  SPLINE  FUNCTIONS  SXCZ)  ANO  SYCZ)  ARE  GIVEN  IN  THEIR  B-SPLINE  10T0 

REPRESENTATION  CB-SPLINE  COEFFICIENTS  CXCJI.RESP.  CYC J) ,J*1 ... .N-K-l)  1080 

ANO  CAN  BE  EVALUATEO  BY  MEANS  OF  FUNCTION  DEKIV.  1090 

CALLING  SEQUENCE:  1100 

CALL  PARAHCX,Y,Z,W,H,Z9,2E,K,S,N,T,CX,CY,FP,I0PT,IPAR,IER)  1110 

1120 

INPUT  parameters:  1130 

X  :  ARRAY. LENGTH  M,  CONTAINING  THE  ABSCISSAE  OF  THE  OATA  POINTS  11*0 

Y  :  ARRAY, LENGTH  M,  CONTAINING  THE  ORDINATES  OF  THE  OATA  POINTS  1190 

W  :  ARRAY. MINIMUM  LENGTH  M, CONTAINING  THE  WEIGHTS  WCI).  1160 

H  :  INTEGER  VALUE. CONTAINING  THE  NUMBER  OF  DATA  POINTS.  1170 

K  :  INTEGER  value, CONTAINING  THE  DEGREE  OF  SXCZ)  AND  SYCZ).  11B0 

S  :  REAL  VALUE. CONTAINING  THE  SMOOTHING  FACTOR.  1190 

IOPT  :  INTEGER  value  WHICH  takes  THE  VALUE  0  OR  1.  1200 

IOPT-O:  the  ROUTINE  WILL  RESTART  ALL  COMPUTATIONS-  1210 

I0PT«1 :  THE  ROUTINE  WILL  START  WITH  THE  KNOTS  FOUND  AT  THE  1220 

LAST  CALL  OF  THE  ROUTINE.  IF  I0PT*1  THE  OUTPUT  1230 

PARAMETERS  T  ANO  N  ARE  INPUT  PARAMETERS  AS  WELL.  12*0 

IF  10P T»1  THE  USER  MUST  PROVIDE  WITH  A  COMMON  BLOCK  1250 

CQMHON/OPT1/NROATACNEST),FPO.FPOLO.NPLUS  1260 

IPAR  :  INTEGER  flag.  1270 

IPAR  *  0:  FOR  EACH  OATA  POINT  CXCD.YCD)  THE  PROGRAM  AUTOMATICALLY  12B0 
CHOOSES  A  CORRESPONDING  VALUE  OF  THE  PARAMETER  Z.  I.E.  1290 

ZC1)«0:ZCI)*ZCI-1)HS0RTCCXCI)-XCI-1))»»2«CYCI)-YU-1)>*P2)  1300 
THE  BOUNDARIES  FOR  THE  PARAMETER  Z  ARE  CHOSEN  AS  FOLLOWS  1310 
ZB  »  ZCI)  :  ZE  *  ZCH).  13Z0 

IPAR  «  t:  THE  USER  HIMSELF  PROVIDES  WITH  THE  VALUES  OF  THE  1330 

PARAMETER  Z  ANO  WITH  THE  BOUNDARIES  ZB  AND  ZE.  13*0 

z  :  ARRAY, LENGTH  M,  CONTAINING  THE  VALUES  OF  THE  PARAHETER  Z  1350 

CIPAR  *  1)  1360 

ZB.ZE:  REAL  VALUES.  CONTAINING  THE  BOUNDARIES  OF  ThE  PARAHETER  Z  1370 

CIPAR  a  1).  1380 

1390 

OUTPUT  PARAMETERS:  1*00 

T  :  ARRAY, LENGTH  NEST  CSEE  DATA  INITIALIZATION  STATEMENT),  1*10 

WHICH  CONTAINS  THE  POSITION  OF  THE  KNOTS. I.E.  THE  POSITION  1*20 

OF  THE  INTERIOR  KNOTS  T C* ♦ 2 ) , . . . T t N-K- 1  )  ,  AS  WELL  AS  THE  1*30 

POSITION  OF  THE  KNOTS  TC 1 ) «T C 2 )* . . . * T C K* 1 ) *  ZB  ANO  ZE  *  1**0 

TCN-K)«...«TCN)  WHICH  ARE  NEEDEO  FOR  THE  B-SPLINE  REPRESENT.  1*50 

CX.CY:  ARRAYS, LENGTH  NEST,  CONTAINING  THE  B-SPLINE  COEFFICIENTS  1*60 

OF  SXCZ).  RESP.  SYCZ).  1*70 

N  :  INTEGER  VALUE • CONT A I  NX NG  THE  TOTAL  NUH3  E  R  OF  KNOTS.  1*30 

FP  :  REAL  VALUE,  WHICH  CONTAINS  THE  S UM C W I*C X I- S X C Z I ) )**2 )  1*90 

*  SUM(WI*<YI-SY(ZI))**2),  i-1, 2, ...M.  1500 

IER  :  ERROR  CODE  1510 

1ER*0:  NORMAL  RETURN.  1520 

ISR»-1:nORMAL  RETURN, SXCZ)  AND  STCZ)  ARE  I NTF RPDL AT ING  SPLINES  1530 
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Source  Listing 


10-0*e-196*  13:53:33 
6-D«e-1984  17:38:13 


WAX-11 

_orao: 


n  :«  g:  h:*o; 

FOR  C  :>0  TO  200  00  BEGIN 
FOR  R  :«0  TO  63  00  BEGIN 

IF  CFCR,C3«25S)  A NO  CN*0)  THEN  BEGIN 
Cl  j*  C:  ri  :*r:  n  :«  n*i: 
end: 

IF  tFCR, 255-0*255)  ANO  CM«0>  THEN  BEGIN 
C2  :*  25S-c:  R2  :*  R:  n  :«m*i; 
end; 

ENO: 

eno: 

WRITELNC 'C0L1«',C1.'R0W1*'.R1.  'CQL2* *  ,  C2,  *R0W2* * ,R2 ) : 

eno: 

(*******a»«*?«*««««*o****»3*********) 

(*  MAIN  PROGRAM  *) 

BEGIN 

WRITELNC 'INPUT  CUT  OR  CONLINE  FILENAME*): 

REAOLN(NAME): 

OPEN  CINFILE. NAME. HISTORY  :«  OLD, 

ACC£SS_METNOO  :*SE3USNTIAL, 

RECORD  LENGTH  J *2 56 , HECORO.T YPE  :*FIX£D): 

OPEN  ( OUTF IlE.'TR.DAT', HISTORY  : »NE W, RcCOR D_LENG TH  :*256, 
RECORD  TYPE  :*FIXE0): 

resetcinfilE): 

REWRITE  COUTrILEi; 

R  :«o: 

WHILE  NOT  EOF  CINFILE)  00 
BEGIN 

RE AO  (INF1LE. IMAGECR3); 

FOR  C  :>  0  TO  255  00 

FCR*1.C*13  :*  IHAGECR.C3; 

R  :•  R*1; 
end: 
first: 

r  :«ri:  c  :«ci:  coir  :*i; 
initial: 

while  CCOC2)  00  BEGIN 

store: 

cmove: 

end: 

FOR  IS«0  TO  COUNT  00 

ACRowcii.coLtna  s -2 55 ; 

FOR  R:«0  TO  63  00  BEGIN 
FOR  C: *0  TO  255  00 

OUTF ILE* CC 3  :«  ACR.C3: 

putcoutfilei; 

eno: 

WRITELNC 'NUMBER* '.COUNT)  : 

closecinfile): 

ENO. 
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1 Q-0#c- 1984  13:53:33  VAX-l 
Source  Listing  6-0*c-198*  17:38:13  _DRA0 

ELSE  IF  ( F0*0 )  ANO  <F7=255)  THEN  BEGIN  C:=>C-15 

store:  r:*r-i:  cdir:*o:  eno 
else  BEGIN  C:«C-i;  coir:*3:  eno: 
eno: 

2:  BEGIN 

IF  (F7*0>  AND  CF6*255)  TNEN  BEGIN  C:*C-lt 
COIR:*3:  ENO 

ELSE  IF  CF6*0)  AND  (F5-255)  THEN  BEGIN  R:*R*i: 

store:  c  :*c-i :  coir:*3:  eno 

ELSE  IF  C  F  5*0 )  ANC  <F4*255;  THEN  6EGIN  R:*R*i: 

C0IR:*2;  ENO 

ELSE  IF  <F4*0)  ANC  <F3*2553  THEN  BEGIN  C:*C*l: 

store:  r:*r«i:  coir:*i:  eno 

ELSE  IF  C F 3*0 )  ANO  (F2*255)  THEN  8EGIN  C:*C*l: 

C0iR:*i;  eno 

ELSE  IF  <F2*0)  ANC  CF1«25S>  THEN  BEGIN  P:«R-i: 

store:  c:*c«i:  coir:*i:  eno 
else  begin  r:*r-i:  coir:*o:  eno: 
eno: 

3:  BEGIN 

IF  (F 1 *0 )  ANO  CF0*255)  THEN  BEGIN  R:*R-i: 

CDiR:*o:  Eno 

ELSE  IF  ( F 0*0 }  ANO  ( F7  =  2  55 )  THEN  BEGIN  C:*C-l: 

store:  r:*r-i;  coir:*o:  eno 

ELSE  IF  t F 7*0 )  ANC  <F6*25S)  THEN  BEGIN  C:«C-l: 
coir:*3:  ENO 

ELSE  IF  (F6*0)  ANO  CF5*25S)  THEN  BEGIN  R:*R*i: 

store:  c:*c-i:  coir:*2:  eno 

ELSE  IF  tF 5*0 i  ANO  (F4*255J  THEN  BEGIN  R:*R«i: 

C0IR:*2:  end 

ELSE  IF  (P4*0)  and  (F3*255)  THEN  BEGIN  C:*C*l: 

store;  r:*r«i:  coir:*2:  eno 
ELSE  BEGIN  c:*C*i;  cdir:*i:  eno: 
bnj: 
eno: 
eno; 


PR0CE0URE  INITIAL: 

BEGIN 

FOR  I  :*0  TO  2S5  DO  BEGIN 

colci 3  :*o : 

RouCl]  :*o: 
eno; 

FOR  j:«0  TO  63  00  BEGIN 
FOR  l!*0  TO  255  00 
ACR,C3:*o: 

eno: 

i:*o: 

eno: 

PROCEDURE  FIRST; 

VAR 

N,M  :INTEG5R; 

BEGIN 
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10-0*c-H84  i3S53:33 
Source  Listing  6-D«e-I9B4  17:38:13 

PROGRAM  TRAC EC  INPUT  » OUTPUT  tlNFILE  > QUTF 1LE) • 

TTPE 

6TTE  «  0..255: 

INAGER0U1  «  PACKED  ARRAY  CO. .2553  OF  BYTE: 

ROU1  -  PACKED  ARRAY  CO. .2573  OF  BYTE: 

VAR 

R,C,F0,Fl,F2,F3,F4,F5.F6,FT,F8  :  BYTE: 

F  :  ARRAY  CO. .653  OF  ROWl: 

A, IMAGE  :  ARRAY  CO. .633  OF  IMAGE  ROU1 : 

INFILE  :  FILE  OF  IMAGEROUI: 

Rl.Cl.R2,C2,CDIK,I, J.COUNT  :  INTEGER: 

RON. COL  :  ARRAY  CO. .5123  OF  INTEGER: 

OUTF1LE  :  file  of  imagerowi: 

NAME  :  PACKED  ARRAY  Cl. .203  OF  CHAR: 

PROCEDURE  STORE: 


BEGIN 

C0LCI3:-C-l: 

R0MCI3:«R-i: 

HRITELNC'COL** ,C0LCI3,*R0u«*,R0MC1 3): 
count  :«i; 
l:«l»i: 
eno: 


PROCEDURE  CMOVE: 


BEGIN 

FO  :«FCR-1,C3 
F4:«FCR*1,C3: 
CASE  COIR  OF 
O:  BEGIN 
IF  CF3-0) 


;  fi:*fcf-i.c*i3:  f2j»fcr»c*i):  F3:«fcr«i.c»i 3: 
FS:»FCR*I,C-13:  F6:»FCR,C-13:  F7:*FCR-1.C-13; 


AND  CF2*255)  THEN  BEGIN  C:*C*I;  CDIR:*l:  ENO 
ELSE  IF  CF2*0)  AND  CF1*255)  THEN  BEGIN  R:»R-i; 

store:  c:*c«i:  coir:*o:  end 

ELSE  IF  CF 1*03  AND  CF0»255)  THEN  BEGIN  R:»R-l; 

coir:«o:  eno 

ELSE  IF  CF0*O)  AND  C  F7*2  5  5  3  THEN  BEGIN  C:*C-i; 

STORE;  r:*r-i;  C0IR:«3:  eno 
ELSE  IF  < F 7*03  AND  (F6*255>  THEN  BEGIN  C: 

CDIR:*3:  end 

ELSE  IF  (F6*C )  ANC  CF5*255>  THEN  BEGIN  R: 

STORE:  c:*c-i :  C0IR:*3;  end 
ELSE  BEGIN  r:*r«i;  c 01 r : *2 :  end; 
eno: 

1:  BEGIN 

IF  (F5«0)  ANO  CF4*255)  THEN  BEGIN  R:»R»i: 


:*c-i ; 


:*R*i: 


COIR : «2 :  END 

ELSE  IF  CF 4* 0 )  AND  CF3«255)  THEN 
STORE!  R:*R+i:  COIR:*l;  ENO 
ELSE  IF  C F 3* 0 3  AND  CF2*2553  THEN 

coir:*i:  eno 

ELSE  IF  ( F2*  0 )  ANC  CF1»255  )  THEN 

store:  c:*c*i:  coir:*o:  end 

ELSE  IF  ( F  1*0)  ANC  CF0*2S5)  THEN 

coir:*o:  eno 


begin  C:»C* 1 ; 

BEGIN  C:*C*i; 
BEGIN  R:*R-i; 
BEGIN  R :* R- 1 ; 


VAK-1 

_ORAO 
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r«*  V  ■VT’.  V.  -r.T. 


I 

\ 

N 

N 

% 

% 

I 


■ 


B 


Source  Listing 


«.  j'c-198*  17:18:40 
4-Dte-U8«  17:18:31 


«  :*  o; 
n  S*  o: 
bright:*  255: 

for  J  :*  0  to  150  do  begin 
for  i :*0  to  63  do  begin 

if <calCit J3*bright)  and  (n*0)  then 
begin 

yl  :*  j; 

*1  :*  i: 
n  :*  n*i; 
end:  (*if«) 

if CcalCi, 255- j3*bright)  and  <e*0)  then 
begin 

y2  :*  255- j; 

«2  :*  i; 

*  :*  *♦  l ; 
end: 

end: C*f or*) 
end  1 C*f oro) 

«ritelnC*«l-',«l,*yl«-,ylt-,2*«,,2,*y2-',v2): 

(*  cut  line*) 

slooe  :*l.o: 
if  *1**2  then 
begin 

slooe  :*  o; 

for  i:**i*i  to  63  do  begin 
for  j:*yl-5  to  y2*5  do 
calCitjl  :*0: 
end ■ <*f oro) 
endlCOife) 
if  slopeOO  then 
begin 

slope  :*<*2-*l J/Cy2-yl); 
if  slope  <  0  then  begin 
slope  J*  abs( slope} ; 
for  j  :*  yi  to  y 2  do  begin 

*  :*  *l*l-roundCCj-yl)bslope>; 
for  i  :»*  to  63  do 

calCi • j3 : *o ; 

end : (Of or*) 
tndi(«ift) 

if  slope  >0  then  begin 

for  j:«  yi  to  y2  do  begin 

*  5*  *l*l»roundCCj-yl)*slope); 
for  i  :**  to  63  do 

calCitJ3  :*0; 
end: (afore) 
end : (Pi f *) 

*riteln( 'step’  *) ; 
endsteife) 

(•put  outfile*) 


for  i :*0  to  63  do  begin 


VAX- 

_0R* 


t  r  f  *unr-  •-  *» 


% 


Sourct  listing 


6-0oc-198<.  17:18:40  VAi- 

6-Doc-1984  17:18:31  „Cn» 


i 


FOR  JS»0  TO  7  DO 

BC64,J  +  2-*8j:*FC63«J-*2483: 
for  i  :*  0  to  63  do  bogin 
fCi+1 • 03  :«  fCi*1.13: 
f Ci  +  1«  2573  :*fti*1.2563: 
ondi 

for  j:«  0  to  255  do  bogin 

fco.  j* 13  :■  f cit J*i 3: 

f 1 65 . j»  1 3  :■  fC64,j*13; 
and: 

fC3.03:*fCi.!3: 
fC0.2573:»fCl,2563: 
f 1 65 , 0 3 : *f C6<  .  1 3 1 
fC65,2573:*U64,2563: 


(•sot  initial  «»i  ain  •) 

■o*:»0; 

•oxi :«o ; 

■in:«255: 

for  i:»  0  to  63  do  bogin 
for  j:«  0  to  255  do  bogin 

dx:*fCi. J3*2#fCi*l. j3*fCi*2, j3-f(i»J*23 
-2*fCi*l, J*23-fCi*2. 3*23; 
dys«fCi*2. J3»2*fCi*2. j*134fCi*2. J*23-fti,J3 
-2*fti»j-»13-fCi«j*23» 
sobalCi.  J3:«round(  (d«**2‘»dy**2)**0.5)  i 
if  >a«  <  »ob%lCi > J 3  than 
u<  :*  $obolCi«43« 
if  oin  >  sob«lCi«J3  than 
■in  :«  sobolCitj3: 

•ndS 

ondl 

rang#  :•  aax~ain; 

<•  rascals*) 

for  l‘>  0  to  63  do  bogin 
for  j:»  0  to  255  do  bogin 

calCi«J3  :*  roundC ( ( sobolCi. J3-*in)*255)/rango) • 
if  <J<«NUH1)  or  ( J>»255-NUM2)  than 
calCi.j)  :»o: 

if  (i<«NUM3)  or  (i>«63-NtjN4)  than 
calCi«j3  :•  0: 

(•  cut  thsshold*) 

IF  THESOO  then  begin 
if  calCiiJ)  <»thos 
than 

colCi»33  i«  0 

•  1  so 

calCit J3  :«  255: 

end: 

and; 

and: 

(•find  point  at  raa  and  aft*) 
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im4  ^ »•  . *  .'-t  /-**  ^ 


OPRINECI)  *  OIAGCI)  *590 

Cx(I )  »  RXCI)  *600 

CYCI)  -  RTCI)  *610 

Gtl. XI)  *  0.  *620 

00  260  J»1.K  *630 

GCI.J)  «  A  (I . J)  *6*0 

260  CONTINUE  *650 

30  300  IT* 1  * N8  *660 

C  THE  ROW  OR  MATRIX  8  IS  ROTATED  INTO  TRIANGLE  BY  GIVENS  TRANSFORMATIONS.  *670 

00  270  1*1. K2  *680 

H( I )  *  BCIT.I)  *690 

2T0  CONTINUE  *700 

XI  •  0.  *710 

Y I  «  0.  *720 

UI  *  PINV  *730 

00  290  J* I T * NK1  *7*0 

IFCWI.EO.O.)  GO  TO  300  *7S0 

PIV  *  MCI)  *760 

C  CALCULATE  THE  PARAMETERS  OP  THE  GIVENS  T R ANSF3RHA TI ON.  *770 

CALL  COSSXNCPIV.WI.DPRIM£CJ).COS.SXN)  *780 

C  TRANSFORMATIONS  TO  RIGHT  HAND  SIDES.  *790 

CALL  ROTATECPIV.CQS.SIM.XI.CX(J))  *800 

CALL  ROTATECPIV.COS.SIN.YI.CYCJ))  *810 

IFC4.EQ.NK1)  GO  TO  300  *620 

12  *  K1  *830 

IFCJ.GT.N8)  12  ■  NK1-J  *8*0 

00  280  1*1.12  *8S0 

C  TRANSFORMATIONS  TO  LEFT  HANO  SIDE.  *860 

CALL  R3TATECPIV,C0S.SIN,MCI*1),6CJ.I))  *870 

MCI)  *  HCI-H  )  *880 

280  CONTINUE  *890 

MC 12*1  )  *  0.  *900 

290  CONTINUE  *910 

300  CONTINUE  -  *920 

C  BACKWARD  SUBSTITUTION  TO  OBTAIN  THE  B-SPLINE  COEFFICIENTS.  *930 

CALL  BACKCG.CX.NKl.K.CX)  *9*0 

CALL  BACKCG.CY.NK1.K.CY)  *950 

C  COMMUTATION  00  PC*}.  *960 

FP  "  0.  *970 

L  «  K2  *980 

00  330  IT* 1 , M  *990 

IFCICin.LT.TCL)  .OR.  L.GT.NK1)  GO  TO  310  $030 

L  *  L«1  5010 

310  LO  *  L-K2  5020 

TERM!  *  0.  5030 

TERM2  *  0.  50*0 

00  320  JM.K1  5050 

LO  *  L 0 ♦  1  5060 

TERMl  *  TERN1*CXCLO)*OCIT. J)  5070 

TERM2  *  TERM2+CTCLO)*OC IT.J)  5090 

320  CONTINUE  5090 

FP  *  FP*WCIT)*CCT£RMI-XCIT))*f2«CTERM2-YCIT))*62)  5100 

330  CONTINUE  JUO 

C  TEST  WHETHER  THE  APPROXIMATION  X* SX PC I > . T* ST P C 2  )  IS  AN  ACCEPTABLE  5120 

C  SOLUTION.  5130 

FPMS  *  FP-S  51*0 

IFCA8SCPPMS).LT.ACC)  GO  TO  **0  5150 

C  TEST  WHETHER  THE  MAXIMAL  NUMBER  OF  ITERATIONS  IS  REACHED.  5160 

IFCITER.EO.MAXIT)  GO  TO  *00  5170 

C  CARRT  OUT  ONE  MORE  STEP  Or  THE  ITERATION  PROCESS.  5190 

P2  *  P  5190 
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F2  *  FFMS 

5200 

-J 

■V] 

IFC ICHECX.NE.O)  GO  TO  3«0 

52X0 

IFCCFZ-F3).GT.ACC)  GO  TO  335 

5220 

‘.Vj 

c  our  initial  choice  of  p  is  too  large. 

5230 

p  *  pao.ie-02 

5240 

t  •  m, 

.VJ 

P3  *  P2 

5250 

F3  *  F 2 

CO  TO  350 

5250 

5270 

335  IFC<F1-F23.GT.ACC)  GO  TO  340 

52  BO 

C  OUR  INITIAL  choice  of  p  IS  too  SMALL 

5290 

..  •  4 

TYPE  4, 'VALUE  OF  P*,P 

P  *  P40.1E*04 

5300 

*•  4 

PI  *  P2 

5310 

FI  »  F2 

5320 

GO  TO  350 

5330 

W **— 4 
► 

C  TEST  WHETHER  THE  ITERATION  PROCESS  PROCEEDS  AS 

theoretically 

5340 

C  EXPECTED. 

5350 

340  IFCF2.GE.F1  .OR.  F2.LE.F3)  GO  TO  4X0 

5350 

ICMECK  *  I 

5370 

C  FIND  THE  NEW  VALUE  FOR  P. 

5330 

P  *  RATION(Pl,Fl,P2,F2tP3tF3) 

5390 

350  CONTINUE 

5400 

C  ERROR  CODES  AND  MESSAGES. 

5410 

u. 

400  IER  *  3 

5420 

GO  TO  440 

5430 

410  IER  *  2 

5*40 

GO  TO  440 

5450 

420  IER  *  1 

5460 

GO  TO  440 

5*70 

430  IER  *  -1 

5480 

\m  * 

440  RETURN 

5*90 

P"" 

ENO 

5500 

( 

SUBROUTINE  BSPUlNU  ,N,X,X,L,H) 

55 10 

C  SUBROUTINE  BSPLIN  EVALUATES  THE  C**l>  NQN-2ER0 

B-SPLINES  OF 

5520 

\V 

C  OEGREE  K  AT  TCL)  <■  X  <  T(l*»  USING  THE  STABLE  RECURRENCE 

5530 

C  RELATION  of  oe  boor  and  COX. 

55*0 

*•  V 

C  THE  DIMENSION  SPECIFICATIONS  OF  THE  FOLLOWING 

ARRAYS  MUST  BE 

5550 

V, 

C  AT  LEAST  MCK*1).MHCK). 

5560 

**.• 

DIMENSION  T(N),HC6).HHC5> 

5570 

ip m 

HU)  «  X. 

00  20  J*1,A 

5580 

5590 

c 

00  XO  1*1, J 

5600 

**  ‘ 

HM Cl)  *  MCI) 

5610 

XO  CONTINUE 

5620 

,**' 

MCI)  *  0. 

5630 

,’t' 

DO  20  1*1, J 

S6-0 

•,* 

LI  »  L  ♦  I 

5650 

LJ  «  LI-J 

F  »  MMCI)/CTCLI)~TCLJ)) 

5660 

5670 

MCI)  ■  HCI)*F*CT(LI)-X) 

5680 

V 

MCIU)  «  F*CX-TCLJ>) 

5690 

,*w 

20  CONTINUE 

5700 

• 

RETURN 

3710 

ENO 

5720 

SUBROUTINE  COSSINCPIV.WI.WW.COS.SIN) 

5730 

C  SUBROUTINE  COSSIN  CALCULATES  THE  PARAMETERS  OF 

A  GIVENS 

57*0 

C  TRANSFORMATION  WITHOUT  SOUARE  ROOTS. 

5750 

STORE  «  PIV4WI 

5760 

00  *  WW*STORE*PIV 

IF(ABSC00;.LT.l.E-36)  DD* 1 . E -36 

COS  ■  WW/OD 

f  9  ?  M 

I 
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f 

t 


SIN  ■  STORE/OO 

UU  *  00 

LI  I  *  C0S*WI 

RETURN 

END 

SUBROUTINE  ROTATEIPIV.COS.SIN.A, B) 

C  SUBROUTINE  ROTATE  APPLIES  A  GIVENS  ROTATION  TO  A  ANC  B. 

STORE  >  8 

B  »  COS*STORE«SIN*A 
A  *  A— PI VPSTORE 
RETURN 
ENO 

SUBROUTINE  BACKIA.Z.N.X.C) 

C  SUBROUTINE  BACK  CALCULATES  THE  SOLUTION  OF  THE  STSTEH  OF 
C  EQUATIONS  ASC  *  l  WITH  A  A  N  X  N  UNIT  UPPER  TRIANGULAR  MATRIX 
C  OF  BANOWlOTH  K«I. 

C  ATTENTION:  THE  FIRST  DIMENSION  SPECIFICATION  OF  MATRIX  A  MUST 
C  BE  THE  SAME  AS  In  THE  CALLING  PROGRAM. 

OIMENSION  AC200,K) .ICNJ.CCN) 

C(N)  *  I(N> 

I  *  N-I 

IFII.EO.O)  GO  TO  30 
00  20  J*2,N 
STORE  «  2C I  ) 

U  •  K 

IF(J.lE.K)  II  *  J-l 
M  ■  I 

00  10  L*1.I1 
M  «  H.  1 

STORE  *  ST0RE-C(M)*ACI.L3 
10  CONTINUE 

CCI)  «  STORE 
I  »  I-l 
20  CONTINUE 
30  RETURN 
ENO 

SUBROUTINE  NKNOT<X,M,T,N,FPI NT, NROATA, NRINT) 

C  SUBROUTINE  NXNOT  LOCATES  AN  ADDITIONAL  KNOT  FOR  A  SPLINE  OF  DEGREE 
C  K  AND  ADJUSTS  THE  CORRESPONDING  P AR AHETER S , I . E. 

C  T  J  THE  POSITION  OF  THE  KNOTS. 

C  N  J  THE  NUMBER  OF  KNOTS. 

C  NRINT  :  THE  NUMBER  OF  KNOTINTERVALS. 

C  F  PI  NT  :  THE  SUM  OF  SQUARES  OF  RESIDUAL  RIGHT  HAND  SIDES 

C  FOR  EACH  KNOT  INTERVAL. 

C  NROATA:  THE  number  OF  OATA  POINTS  INSIDE  EACH  KNOT  INTERVAL. 

C  THE  ARRATS  T.FPINT  AND  NROATA  MUST  HAVE  THE  SAME  DIMENSION 
C  SPECIFICATIONS  AS  IN  THE  CALLING  PROGRAM. 

OIMENSION  X<H),T(200),FPINT<200),NROATAC20Q) 

K  ■  (N-NRINT-1 )/2 

C  SEARCH  FOR  KNOT  INTERVAL  T ( NUMBER** )  <■  X  <■  T t NUMB E R *K ♦ 1 )  MMERE 
C  FPlNT(NUMBER)  IS  MAXIMAL  ON  THE  COnOITION  THAT  NRDATACNUMBER) 

C  NOT  EQUALS  zero. 

FPMAX  »  0. 

JBEGIN  *  1 
00  20  J*l, NRINT 

JPQINT  «  NROATA(J) 

IFCFPNAX.GE.FPINTC J)  .OR.  JPOINT.EO.O)  GO  TO  10 
FPMAX  »  FPINTCJ) 

NUMBER  *  J 
MAXPT  »  JP  01  NT 
MAXBEG  *  JBEGIN 


10  JBEGIN  »  JBEGIN+JPOINT  *1  6400 

20  CONTINUE  6410 

C  LET  COINCIDE  THE  NEW  KNOT  T ( NUN  BE R »K *1 )  WITH  A  DATA  POINT  X(NRX>  6420 

C  INSIDE  THE  OLO  KNOT  INTERVAL  T<NUN8£R*K>  <*  X  <«  T ( NUMB ER+K ♦ 1 ) .  6430 

IHALF  *  MAXPT/2«1  6440 

NRX  *  MA  X  BE  G* I HALF  64  SO 

NEXT  x  NUH8ER«1  6460 

IF(NeXT.GT.NRINT)  GO  TO  40  6470 

C  ADJUSTS  THE  DIFFERENT  PARAMETERS.  6460 

00  30  JxNE XT . NRI NT  6490 

JJ  *  NEXT+NRINT-J  6500 

FPINT(JJ*1)  *  FPINUJJ)  6510 

NRDATA(JJ*1)  «  NROATA(JJ)  6520 

JK  x  JJ*K  6530 

T(JK*1>  «  T  (  JK  )  6540 

30  CONTINUE  6550 

40  NRDATA(NUHBER)  «  IHALF-1  6560 

NRDAT  A (NE  XT )  «  HAXPT-IHALF  6570 

FP INT {  NUMBER )  *  FPMA X*FL DA T< NR D AT A < NUHB E R ) >/PLDA T ( HA  X P T )  6590 

FPINT(NEXT)  *  FPMAX*FLOAT(NRDATA(NEXT>  3/FLOATCMAXPT)  6590 

JK  *  NEXT«K  6600 

TCJK)  *  X{ NR  X }  6610 

N  *  N»1  6620 

NRINT  «  NRINT* 1  6630 

RETURN  6640 

ENO  6650 

SUBROUTINE  OISCOCT ,N, K2, B>  4660 

C  SUBROUTINE  DISCO  CALCULATES  THE  DISCONTINUITY  JUMPS  Dc  THE  KTH  6670 

C  DERIVATIVE  OF  THE  B-SPLINES  OF  DEGREE  K  AT  THE  KNOTS  T ( K ♦ 2) . . T C N-K- 1 >  6690 

C  THE  FIRST  DIMENSION  S°ECI FI CATI ON  OF  THE  MATRIX  E  MUST  BE  THE  SAME  AS  6690 

C  IN  THE  CALLING  PROGRAM:  M  MUST  HAVE  A  DIMENSION  SPECIFICATION  AT  6700 

C  LEAST  2*K*2.  6710 

DIMENSION  TCN),BC200.K2).HU23  6720 

K 1  x  K2-1  6730 

K  *  Kl-1  6740 

NKI  *  N-Kl  6750 

00  40  L*K2,NK1  6760 

LMK  *  L-Kl  6770 

30  10  J=l.Kl  6780 

IK  «  J«KI  6790 

LJ  «  L*J  6800 

LK  *  LJ-K2  6810 

HCJ)  =  T(L)-T(LKJ  6820 

HCIK)  X  T(L3-T(LJJ  6830 

10  CONTINUE  6840 

LP  x  LMK  6350 

00  30  Jx 1 1  K2  6860 

JK  x  j»K  6870 

PROO  x  1.  6880 

00  20  IxJ.JK  6890 

PROO  «  PR004MCI)  6900 

20  CONTINUE  o  9 1 0 

LK  x  LP«K1  6920 

BCLMK.J)  x  (T(LK)-T(LP) J/PROO  6930 

LP  «  LP*1  6940 

30  CONTINUE  6950 

40  CONTINUE  6960 

RETURN  6970 

END  6980 

FUNCTION  RATION (PI *Fl,P2tF2fP3*F3)  6990 

C  GIVEN  THREE  POINTS  (P1»F1),(P2«F2)  ANO  (P3.F3),  FUNCTION  RATION  7010 
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C  GIVES  THE  VALUE  OF  P  SUCH  THAT  THE  RATIONAL  INTERPOLATING  FUNCTION  7010 

C  OF  THE  PORH  RCP  )  «  CU*P*V)/CP*W>  EQUALS  ZERO  AT  P.  7020 

IFCP3.GT .0. )  GO  TO  10  7030 

C  VALUE  OF  P  IN  CASE  P3  »  INFINITY.  70*0 

P  «  CPl*CFl-F3)*F2-P2*(F2-F3)pFl>/aFl-F2)*F3)  7050 

GO  TO  20  T060 

C  VALUE  OF  P  IN  case  P3  "■  INFINITY.  7070 

10  HI  *  F1PCF2-F3)  7080 

H2  ■  F2*CF  3-FI )  7090 

H3  ■  F3*(F1-F2>  7100 

p  «  -CP1»P2*H3*P2*P3*H1*P3PP1»H2)/<P1*H1«P2PH2*P3PH3>  7110 

C  AO JUST  THE  VALUE  OF  P1.F1.P3  ANO  FJ  SUCH  THAT  FI  >  0  AND  F3  <  0.  7120 

20  IFCF2.LT.0.)  GO  TO  30  7130 

PI  *  P2  71*0 

FI  »  F 2  7150 

GO  TO  *0  7150 

30  P3  ■  P2  7170 

F3  *  F2  7180 

*0  RATION  *  P  7110 

RETURN  7200 

End  7210 

FUNCTION  OERIVO.N.C.NKl.NU.ARG.L)  7220 

C  FUNCTION  OERIV  EVALUATES  A  SPLINE  SCX)  OF  DEGREE  K  WHICH  IS  7230 

C  GIVEN  IN  ITS  NORHAL 1 ZED  8-SPLINE  REPRESENTATION  DR  CALCULATES  72*0 

C  DERIVATIVES  OF  ANT  SPECIFIED  ORDER  NU.  7250 

C  7250 

C  CALLING  SEQUENCE  7270 

C  VALUE  *  OERIVCT.N.C.NM.NU.ARG.L)  7230 

C  7290 

C  INPUT  PARAMETERS:  7300 

C  T  :  ARRAY, MINIMUM  LENGTH  N,  WHICH  CONTAINS  THE  POSITION  7310 

C  OF  THE  KNOTS  OF  SCX),I.E.  THE  POSITION  O'  THE  INTERIOR  7320 

C  KNOTS  T(K*2),...TCN-K-13  AS  WELL  AS  THE  POSITION  OF  THE  7330 

C  KNOTS  T(1),...T(K«1)  ANO  T CN-* ) , . . .T CN  )  WHICH  ARE  NEEDED  73*0 

C  FOR  THE  0-SPLINE  REPRESENTATION.  7350 

C  N  :  INTEGER  VALUc  GIVING  THE  TOTAL  NUMBER  OF  KNOTS  OF  SIX).  7350 

C  C  :  ARRAY, LENGTH  MAI,  CONTAINING  THE  0-SPLINE  COEFFICIENTS.  T370 

C  NK1  :  INTEGER  VALUE. GIVING  THE  OIMENSION  OF  S(X),I.E.NK1  >  N-K-l.  T380 

C  NU  :  INTEGER  VALUE  WHICH  GIVES  THE  OROER  OF  THE  DERIVATIVE.  T390 

C  ARG  :  REAL  VALUE, GIVING  THE  VALUE  OF  THE  ARGUMENT.  T*00 

C  L  :  INTEGER  VALUE  WHICH  SPECIFIES  THE  POSITION  OF  THE  ARGUMENT  7*10 

C  X.E.  TCL)  <*  ARG  <  TCL+l)  OR  7*20 

C  L  *  NX1  IF  ARG  *  T(NK1«1).  7*30 

C  7**0 

C  OUTPUT  PARAMETER:  7*50 

C  VALUE:  REAL  VALUE, GIVING  THE  VALUE  OF  THE  NUTH  DERIVATIVE  OF  7*60 

C  SCX)  AT  X  *  ARG.  7*70 

C  7*80 

C  OTHER  SUBROUTINES  REQUIRED:  NONE.  7*90 

C  7500 

C  RESTRICTIONS:  7510 

C  1)  NU  >*  0  7520 

C  2)  TCK«1)  <«  ARG  <»  TCNM  ♦ 1 )  T530 

C  THE  OIMENSION  SPECIFICATION  OF  THE  ARRAY  H  MUST  8E  AT  LEAST  K*l.  75*0 

OIMENSION  TCN),C(NK1),HC6)  7550 

OERIV  *  0.  7560 

K 1  «  N-NK1  T5  TO 

IFCNU.LT. 0  .OR.  NU.GE.K1)  RETURN  7580 

00  100  I  *  1 , K 1  7590 

IK  *  L ♦! -K I  7600 

H  C I )  »  CCIK)  ’*  ’  P 
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.v  <.• '.r-rv 


100 


200 

300 


400 

500 


600 


CONTINUE 

lECNU.EC.O)  GO  TO  300 
NU1  *  NU*1 
00  200  J*2,NU1 
00  200  JJ*J,K1 
I  «  J*H1-JJ 
Cl  «  LM-K1 
LJ  *  C*I-J»1 

H<1)  ■  <h(I)-h(I-1 ))/<T<LJ) 

continue 

IE<NU.EC.K1-1)  go  to  500 
NU2  «  NU*2 
OC  400  J *NU2  •  K  1 
00  400 

I  «  J««1-JJ 
Cl  «  C«l-*1 
CJ  «  C*I-J*1 


-TClXJ) 


Htl) 

continue 


(U«C*T(U))»H(I).(TCU)-HG)»K(M))/(T(U,.ULI)) 


OERIV  *  H(ll) 

IFCNU.EC.O)  RETURN 
00  600  1*1. NU 

DERI*  *  OERIV*ELO*T<H-I) 
CONTINUE 
RETURN 
ENO 


T  6  2  0 

76  30 

7640 

7650 

7660 

7670 

7680 

7690 

7700 

7710 

7720 

7730 

7740 

7750 

7760 

7770 

7780 

7  7  »0 

7800 

7810 

7820 

7830 

7860 

7850 

7860 

7870 


•s 
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Sourcs  Listing 


10-D«c-U34  16:37:51 
10-D*c-19B4  16:37:39 


VA  1-1 
_DRA0 


PROGRAM  BSPL I NEC  INPUT, OUTPUT, I NFILE.OUTFILE) ! 


TYPE 

BYTE  *  0..2SS: 

I  HAGER  DW 1  *  PACKED  ARRAY  C0..25S3  OF  BYTE: 

ROU1  ■  PACKED  ARRAY  CO. .2573  OF  BYTE! 

SHIP  »  PACKED  ARRAY  C 0 . . 63 , 0. . 25 5 3  OF  BYTE; 

DA 1  *  PACKED  ARRAY  CO. .5123  OF  REAL! 

DA2  *  PACKED  ARRAY  Cl. .3003  OF  REAL! 

VAR 

R. C,F0,F1,F2,F3,F4,F5.F6.F7,F8  :  BYTE; 

F  :  ARRAY  CO. .653  OF  ROWi; 

A, IMAGE  :  SHIP; 

INFlLE  ;  file  OF  IMAGEROUl: 

SPX,SPY,SPX1,SPY1  ;  PACKED  ARRAY  CO. .5123  OF  REAL! 
R1,C1,R2,C2,CDIR,I, J ,J1 . COUNT , NEG  : INTEGER: 
tempx.tempy.ra  :  real: 

N,I0PT,K.IPAR,N,IER,ANS,NU,ANS1,ANS2.ANS3  :  integer: 
NK1.NEN0.L.MET  :integer; 
y ,  Z  :PACkEO  ARRAY  C0..5123  of  real: 

S. ZB.ZE.FP  :  real: 

ARG,THETA,T0LE  :  REAL: 

FLAG_BEGIN,plaG_ENO.AK,LUMP  i  integer; 
beg, In  :  packed'array  ci.,53  of  real: 

SEGN.ENN  :  PACKED  ARRAY  Cl. .53  OF  INTEGER: 

T. CX.CY  : PACK  ED  ARRAY  Cl. .3003  OF  REAL : 

COL, ROW  : P ACKE  D  ARRAY  CO. .5123  OF  REAL! 

X,Y  :PACKEO  ARRAY  C0..5123  OF  real: 
OCY.OCX.ARE.CYl.DMAX.DMIN  :  REAL; 

CVMIN.CYMAX.HAxCY  :  PACKEO  ARRAY  Cl. .53  OF  REAL.* 
AREA  :  PACKEO  ARRAY  Cl. .53  OF  REAL: 

outfile  :  file  of  imagerohi: 

NAME  :  PACKEO  ARRAY  Cl. .203  OF  CHAR; 

PEAK, STAR. TER  : PACKED  ARRAY  Cl. .S3  OF  REAL: 


c*  filter  the  points  *> 

PROCEDURE  STORE; 

BEGIN 

tempx:*c-i:  tempy:*6«.-CR-d: 

IF  I>0  THEN  BEGIN 

FOR  j:*0  TO  H  00  BEGIN 

IFCCOLCJ J*TEHPX)  AND  < ROW C J 3 * TE HP Y ) 
THEN  I :> j; 

end: 

C0LCI3  :*TEmpx; 
rowcI3  :*tempy: 
m:«i  : 

I :  «I *1 : 

•no 

ELSE  BEGIN 

COlc 1 3  :«tempx: 

Rowcn  :*tempy: 
m:»o; 
i  :«i: 
eno: 
end: 


!_ "  .  •  *  '  a  '  •  *  .  *  »  .  *  .  '  i  *it  **  •  "  -  ' 
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PROCEOURE  cnove: 

BEGIN 

FO  :»FCR-1,C3:  F1:*FCR-1.C*13:  F2  :*FER ,C *1 3 S 
F3:«FCR*l,C«n;  F4:*F£R*I  ,C3!  F5  S  »F£  »♦  1  ,  C-l  3  : 
F6:*FtR»c-n:  F7:«fcr-i,c-i3: 

CASE  COIR  OF 
o:  BEGIN 

IF  ( F  3*0 )  ANO  (F2*255>  THEN  BEGIN  C:*C«1; 

COI R : *1 ; 

END 

ELSE  IF  (F2*0)  ANO  (Fl*255)  THEN  BEGIN  R:*R-i: 

store:  c:*c*i:  coir:*o:  eno 

ELSE  IF  (FI *0 1  ANO  (F0*255)  THEN  BEGIN  R:*R-l; 
COIR:*o ;  END 

ELSE  IF  (F0»0)  AND  (F7*255)  THEN  BEGIN  C:*C-i: 

store:  r:*r-i;  coir:*3:  end 

ELSE  IF  (F7*0)  And  (F6*255)  THEN  BEGIN  C:*C-l: 
COI R : • 3 ;  END 

ELSE  IF  (F6*0)  ANO  (F5*255)  THEN  BEGIN  R:*R«i; 

store;  c:«c-i :  coir:*3:  end 

ELSE  BEGIN  R:*R«l;  COIR:*2:  END! 

eno: 

l:  BEGIN 

IF  ( F 5*0 J  ANO  ( F4*255 }  THEN  BEGIN  R:*R»i: 
C0IR:*2;  ENO 

ELSE  IF  (F**0)  ANO  (F3*2553  THEN  BEGIN  C:*C*1 

store:  r:*r«i:  coir:*i:  eno 

ELSE  IF  ( F 3*0 1  AnC  (F2*255>  THEN  BEGIN  ti*C*l 
COIR: *1 ;  ENO 

ELSE  IF  (F2*0)  ANO  <F1*25S)  THEN  BEGIN  R:*R-1 

store;  c:*c«i:  cdir:*o;  eno 

ELSE  IF  (FI *0 )  ANO  (F0-255J  THEN  BEGIN  R:*R-l 
CDIR:*0;  ENO 

ELSE  IF  ( F0*0 3  ANO  (F7*2S5)  THEN  BEGIN  C:*C-1 

store:  r:*r-i;  coir:«o;  eno 
else  begin  c:»c-i;  coir:*3;  eno; 
eno: 

2:  BEGIN 


•  IF  (F7*0>  ANO  ( F 6* 2 5 5 )  THEN  BEGIN  C:«C-l; 
cdir:*3:  Eno 

ELSE  IF  (F6*0)  ANO  (F5*255)  THEN  BEGIN  R:*R*i; 

store;  C:*c-i;  coir:*3;  eno 

ELSE  IF  ( F 5*0 1  ANO  (F4*255)  THEN  BEGIN  R:*R*l; 
C0IR:*2 :  ENO 


ELSE  IF  { F4*0 >  ANO  (F3»25S)  THEN  BEGIN  C:*C*i: 

store:  r:*r*i;  coir:*i;  end 

ELSE  IF  (F3*0)  And  (F2*255l  THEN  BEGIN  C:*C«i: 
C0ir:*i;  eno 

ELSE  IF  (F2*0 1  AND  (Fl*255)  THEN  BEGIN  R:*R-i; 

store:  c:*c«i:  coir:*i;  eno 

ELSE  BEGIN  R :  *R-1 ;  C0IR:*0;  ENO: 

end: 

3:  BEGIN 

IF  ( F 1 *0  )  AND  (FO*  255  )  THEN  BEGIN  R:»R-i: 
C0IR:*0;  ENO 

ELSE  IF  ( F  0  *  0 )  ANC  (F7*255)  THEN  BEGIN  C:«C-l: 
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VA  «-  1 
_DRA0 


store:  r:»r-i;  coir:*o:  end 

ELSE  IB  CF7-0)  ANO  <F6=255)  THEN  BEGIN  C:*C-i: 
C0IR:»3:  eno 

ELSE  IF  CF6»0>  ANO  CF5*255)  THEN  BEGIN  R:*R*l: 

store:  c:«C-i:  cdir:»2:  eno 

ELSE  IF  (F5«0)  ANO  (F**255>  THEN  BEGIN  R:*R*i: 

coir:>2:  eno 

ELSE  IF  (  F<»*0  )  ANO  (P3*25S>  THEN  BEGIN  C:*C*l: 

store:  r:«r*i:  coir:«2:  eno 

ELSE  BEGIN  c:«C«i;  cdir:*i:  eno: 

eno: 

eno: 

eno: 


proceoure  initial: 

BEGIN 

FOR  I  : • 0  TO  2SS  DO  BEGIN 
XCI3  :*o: 

Ten  :»o; 

eno: 

I  :>0: 
eno: 


PROCEDURE  first: 

VAR 

n.h  :integer; 
begin 

n  :•  o:  h:»o: 

FOR  C  :*0  TO  200  00  BEGIN 
FOR  R  :«0  TO  63  00  BEGIN 

IF  (FCR,C3*2S5>  ANO  CN»0)  THEN  BEGIN 
Cl  :*  c:  ri  :«R;  n  :»  n*i: 

(NO: 

IF  <FCR,255-C3»2S5)  ANO  <H«0)  THEN  BEGIN 
C2  :*  255-c;  R2  :«  r:  n  :«h+i; 
eno: 
eno: 
eno; 
eno: 

procedure  rotate: 

BEGIN 

neg:«o; 

IF  R1«R2  then  begin  theta;«o.o: 

FOR  l:«0  TO  H  00  BEGIN 

XCI3:<COLCI3:  re  I  3 : «RDNC I  3 : 

eno: 

ENO 

ELSE  THETA  :*  ABSC  ARCTANUR2-R  1  )/<C2-Cl  )))  ; 

IF  (THCTAOO)  ANO  <R2>R1)  THEN  BEGIN 
FOR  l:«0  TO  N  00  BEGIN 

XCI3:«C0LCI3*C0S<THETA)-R0MCI3*SIN(THETA); 

TCI3r»COLCI3*SlM(THETA)*ROVCI3*COS(THETA): 
IF  TCI3>>63.0  THEN  neg:>neg*i: 
eno: 
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Vi  X-l 

_orao 


END 

ELSE  IF  ( THE  TAOO  )  AND  (R2<R1)  THEN  BEGIN 
FOR  l:>0  TO  M  00  BEGIN 

XCI J:*COLCI30COS{THETA)*ROWCI30SINCTHETA)  ; 
TCI3:»-COLCI3oSIN(THETA>*ROWEn*COS(TMETA): 

IF  YCI3>»63.0  THEN  NEC : "NEG+l ; 
end: 
end: 

(*  FILTER  *) 
l:*o; 

FOR  KOO  TO  M  00  BEGIN 

tehpxoxcki;  tehpy:«ycx3: 
if  i>o  then  begin 

FOR  j:*0  TO  H  00  BEGIN 

IF  (XCJ3*T£HPX)  and  <YCJ3*TEHPY>  then  ::=j; 
end: 

xci3:*tehpx;  tcI3:«tehpt:  m:»i;  i:=i*i: 

ENO 

else  begin  xcid:*tehpx:  yci3:*tehpy;  h:*o:  i;*i: 
eno: 
eno: 
eno: 

proceoure  parahc  x:oai:  y:oai:  var  z:oai;  u:oai: 
h:integer:  var  zb:real;  var  ze:real;  k:integer; 
S:real:  var  n:integer;  var  t:oa2:  var  cx:oa2: 
var  ct:oa2:  var  ppsreal;  iopt:integer: 
ipar :i nteger :  var  ier:integer);  fortran; 

PROCEDURE  INITTCSPEEO: INTEGER!  :  FORTRAN; 

proceoure  vminoocxhin:real:  xrangesreal;  yhinsreal; 

yrange :real)  :  fortran: 
proceoure  hoveacxjreal;  t:realj:  fortran; 
procedure  orauacxjreal;  t:read;  fortran; 

PROCEDURE  ANCHO(ICHAR:lNTEGER);  fortran; 

PROCEOURE  FINITTCH5INTEGER;  12 :  INTEGER)  :  fortran; 
proceoure  dashacx:real;  t:real;  l:integer):  fortran; 

FUNCTION  OERIVCtREF  T;0A2;  NSINTEGER:  Cx:0A2: 
NKIZINTEGER:  NU :INTEGER:  ARG:REAL; 

l:integer)  :real;  Fortran; 

PROCEOURE  SPLCOEF; 

BEGIN 

i:*s:  luhp:*o; 

WHILE  (I<*N-S)  ANO  <LUHP«0)  00  BEGIN 

iF(crci-n>crci3)  ano  ccyci*i3>cyci3>  ano 
CCYCI*23>CYCI«13)  THEN 

LUHP:«1 ; 

is«i»i: 

eno; 

IF  LUHP* 1  THEN  BEGIN 

Flag.begin:«0;  ax:*i;  fcag_Eno:*q;  i:«5: 
while  ion-5  DO  BEGIN 

IF  FLAG_BEGIN*0  then  begin 

IF  CCYCI-13KTCI3)  ANO  C  CYC  I  ♦  1  3  >C  Y  C I  3  )  ANO 
(CTEI*23>CYCI*13)  THEN  BEGIN 
3EGCAX1:«TCI3;  FL  Ai.SEGIN : » 1 ;  3EGNCAR30I: 
ENO; 
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END 

ELSE  IF  FLAG_B€GIN»1  THEN  BEGIN 

IF  <CT£I-13>*CYCI3>  AND  < C Y C I 3>»C YC I ♦ 1 3 1  AND 
<CTCI*23>«CYCI*13>  THEN  BEGIN 
ENCAX3:»TCI+13:  FLAC_ENDi«l;  ENNC AK 3 {  =  I *1 : 

end: 

end; 

IF  (FLAG_BEGIN*1 )  AND  (FLAG_EN0«1)  then  begin 
flag_BEGIn:*o:  flag_Eno:*o:  ax:«ax*i:  i :  =  i - i ; 
eno: 

I :*l ♦ x  s 

eno: 

END 

ELSE  IF  LUMP»0  THEN  BEGIN 

flag_begin:xo;  flag_end:»o;  i:*5:  ax:=i; 

WHILE  I  <  *N-5  DO  BEGIN 

IF  FLAG.BEGIN*0  THEN  BEGIN 

IF(CYCI-13>*CYCI3>  AND  CCY C I ♦ 1 3 >C YE  I  3 >  AND 
(CYCI*23>CTCI3)  THEN  BEGIN 
3EGCAX3:«TCI3:  FLAG.BE GI N : « 1 :  BEGNE  AX  3 :  *1 ; 

eno; 

ENO 

ELSE  IF  FLAG.BEGINxl  THEN  BEGIN 

IFCCYCI— 13>*CYCI3)  AND  C < C Y E 1 3< *C YE  I ♦ 1 3* T Ol E  1 
AND  CCYCI3>*CYEI*13-TOLE)>  ANO 
(<CYEI3<=CYEI»23*TOLE)  ANO 
CCTCI3>xCTCI*23*TOLE))  THEN  BEGIN 
encak3:xtci3 :  flaG_Enojxi;  enneax3:*i; 
eno: 
eno; 

IF  <FLAG_BEGIN*1>  AND  IFLAG.ENDxi)  THEN  BEGIN 
flag_6EGIn:xo;  flag_Eno:xo;  ak:xak»i; 
eno; 
i:*I*i: 
eno; 

IF  FLAG_BEGINxl  then  begin 
I:«begncax3: 

WHILE  C I <*N-5 )  and  CFLAG_EN0*0)  00  BEGIN 
IF  (CYEI-13>»CYEI3)  ANO 
((CTCI3<xCTCI*13*TOLE) 

ANO  tCTCI3>«CTCI*13-TQLEJ)  THEN  BEGIN 
ENCAX3:«TCI3:  ENNCAK3:»I:  FLAG  END : x 1 ; 

eno; 

l:*l*s: 

end; 

end; 

eno; 

eno; 

PROCEDURE  ALUHP; 

BEGIN 

AX : xi ; 

WHILE  C  E  NNC  Ax  3 ) >0  DO  BEGIN 

areacax3:*o;  cyhaxeax3:«o.o:  cyhineax3:*iooo.o; 
FOR  I  :  * ( 3  EGNC AX3  3  TO  ( ENNC  AX  3  )  00  BEGIN 
IF  CYMAXCAX3  <«  CTCI3  THEN  BEGIN 
CTNAXIAX3SxCTCI3;  HAXCYCAX3:«TEI3 
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eno  : 

IF  CYMINCAX3  >-  C  T  Cl  3  THEN 
ctnincah3:*ctc13: 

eno: 

CTl:*CYCBEGNCAK33-CTMlNlAX3: 

FOR  I  :*<BE&NCAK3)  TO  ( ENNC  AN  3-1 3  00  BEGIN 
0CT:«CTCI«13-CTCI3 ; 

0CX:*CXCI*13-CXCI3: 

AREACAN3:*AREACAN3«CT1«OCX-»OCTSOCX/2.0: 
WRITELNC *AREA»*,AREACAN3.*  I*'. It*  AX*', AO: 

cti:*cti*oct: 

Eno; 

ax:>aui  : 
eno: 
eno: 

PROCEOURE  PCSINX: 

BEGIN 

i:*0;  het:*o: 

WHILE  CI<»H-13  AND  ( ME  T  «0  3  00  BEGIN 
IF  TeHPX=2CI3  THEN  BEGIN 
TEhpt :=xCI3;  mEt:*i: 
eno: 
i  s  *  i  ♦  i : 
eno; 
eno; 

(««»*««»  «*«««««  e*  ««»««*«  *e*s*e  *«**«) 

(*  HA  I N  PROGRAM  «) 

BEGIN 

WRITELNC 'INPUT  COT  OR  CONLINE  FILEnAHE*): 
REAOLNtNAME  3  ; 

OPEN  (INFXLE. NAME. HISTORY  !<  0L0, 

ACCESS  METHOD  :»SEQUENTIAL, 

RECORC.LENGTH  :«256.R£CORO_TTPE  :«FIXE03: 
RESETC INFILE3: 

fl  :«o: 

WHILE  NOT  EOF  (INFILE3  00 
BEGIN 

READ  (INFILE. IHAGECR33; 

FOR  C  :=  0  TO  255  00 

FCR»1.C»13  :»  IHAGECR.C3: 

r  :*  r*i; 

eno: 

close(in*ilE3; 

first: 

p  :«ri:  c  :«ci:  coir  :*i: 
initial; 

while  (COC23  00  BEGIN 
store: 
chove: 
eno; 
rotate: 

FOR  I : *0  TO  H  00  BEGIN 

IF  NEGOO  THEN  T C 1 3  : *  TC I  3-2  0. 0 ; 
eno: 

m:*m*i;  x:*2:  s:*0.i:  iopt:«o;  ipar:»o: 
for  i:*o  T 3  M  00 
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WCI3:=1.0: 

ANS:*1  ; 

WHILE  AN  S3 1  00  BEGIN 

WRITELNCOLO  S*'.S,*NEW  S«')I 
REAOCS): 

WRITELNCOLO  K**, A, •»«'): 

REAO(A)  : 

PARA!i(X,»,Z,W,K,ZB.ZS,X,S,N,T,CX,CY,FP,IOPT, 
IPAR ,IER)  ; 

WRITELNC*  S='.SC  IER=',IER,'  N=',M,'  N=',N. 

*  CXCN3* '.CXCNJ.  '  CTC  N— 4  3=*.CYCN— 4j); 
WRITELNC  PRINTING  YES=1  PLOT  X -Y  =  2  .  C  X  ,  C  Y  =  3 
*X -Y-C  X-C  Y  *4 ' , *  NO  =  5'3: 

RE  AD (ANSI): 

IF  A  NS  1  *  1  THEN  BEGIN 
FOR  I : = 1  TO  N  DO 

WRITELNCCX.'.CXCI).'  CT*',CTCn, 

Te'.Tcn); 

ENO 

else  if  ansi=2  then  begin 

IMTT(960): 

VwlNO0(0.O.2C P-13,0. 0.256.0): 
M0VEAC2t03,XC03)  : 

FOR  I :*0  TO  M-l  00 

DRAuAccn  .xci )) ; 

MOYEAczcai.rcoi): 

FOR  I : *0  TO  M-l  00 

oasma(zci3.t:i3,2); 

finittco.o); 

ENO 

ELSE  IF  ANS1*4  THEN  BEGIN 
INITTO60)  : 

VWINOOCO.O.IIM- 13, 0.0, 256.0): 
MOVEA(ZC03,XI03): 

FOR  I : *0  TO  M-l  00 
0RAWA(ZCI3,XCI3): 

MOVEA(ZC03.TC03)  : 

FOR  I : =0  TO  M-l  00 
ORAWACZCI3.TCI3): 

M0VEA(TC4],CXC4)); 

FOR  l:=4  TO  N— 4  00 

OASMA(TCI3.CXCI3,2): 

FOR  I : « 4  TO  N  — 4  00  BEGIN 

M0VEA(TCn-1.5.CXCI3-l.S): 

ANCHOC  III  )  ; 

enc; 

M0VEA(TC43,CU43): 

FOR  l:»4  TO  N— 4  00 
OASMA(TCI3,CTCI3.2): 

FOR  I : « 4  TO  N-4  oo  BEGIN 
MOVE  A(TCI3-1.5tCTtl)-1.5); 

ANChO(HI)  ; 

ENO; 

Dmin:*o.o:  oma x : *2 56. o :  tempx:=omin; 

WHILE  TEMPX<*ZZM-13  OC  BEGIN 
MOVEA<  TEMPX.OMIN): 

ORAWAC TEMPX.OMIN) ; 


VAX-1 

_qsac 
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Sourc#  Listing 


t 0-0*c- 198  4  16:37:51 
10-D»c-1984  16:37:39 


WAX-1 

_prac 


QRAWACTEMPX .OH  AX )  1 
TENPX:»TEHPX*10.QS 

eno: 

omins-o.os  omax:*icm-13;  tempx:*omin; 
WHILE  TENPXX-256.0  00  BEGIN 

moveacdmin.tempx); 

ORAMACDMIN.TEMPX) S 
QRAWA(DMAX,TEHPX)S 
TEMPX:»TEHPX*10.0S 
ENOS 

FINITTCO.O); 

ENO 

ELSE  IP  ANS1*3  THEN  BEGIN 
INITT (960 ) S 

VWINOOCO.O, TEN), 0. 0.256. 0)S 
M0VEACTC43.CXC43) S 
FOR  15*4  TO  N-4  00 
0RAWA<TCI3,CXCI3)S 
FOR  I  S»4  TO  N-4  00  BEGIN 
M0VEACTCI3-1.5.CXCI3-1.S); 

ANCHOC 1 11 ) ; 

ENOS 

H0VEACTC43.CTC43); 

FOR  1 S»4  TO  N-4  00 
0ASHA<TCI3,CVCI3,2>S 
FOR  I:*4  TO  N-4  00  BEGIN 
MOWEACTCI3-1.5.CTCI3-1.53S 
ANCHOUIDS 
ENOS 

FINITT(O.O): 

eno: 

(*  IMPORTANT  FORTRAN  OECLEAR  FROM  1  *) 

NR  1 1  «N— R- 1 5  LS*R*1S  NENO:»N-RS  J1:«0S 
FOR  I :*L  TO  NENO  00  BEGIN 
ARG:«TCI3S 

WHILSC ARG>*TCL*13)  ANO  CL<NX1)  00 
LS*L*1 S 

SPX1CJ13:*  DERIVCT.N.CX.NRl.O, ARG.L); 
SPUCJ13SX  OERIV  CT  «N.CT  .  NR  1  .O.ARG.L)  S 
Ji:«Jl«ls 

ENO; 

j:*Os  L :*R *1 s 

FOR  i:>L  TO  NENO  CO  BEGIN 
ARGS*TCI3S 

WHILE  <ARC>«TCL«133  ANO  CL<NK13  DO 

l:»l*i: 

tempx:*tci»13-args 

IF  TEHPX>*7.0  THEN  BEGIN 
WHILE  ARG<TCI ♦ 1 3  00  BEGIN 

spxcj3  :«  oerivct.m.cx.nri.o, arg,l>: 

SPTCJ3  S*  OERIVCT.N.CT.NRl, 0.ARG.L3! 

arg: *arg*2.0 ;  j:>j»i; 
eno; 

ENO 

ELSE  BEGIN 

SPXCJ3  S«  OERIVIT.n.CX.NRl, O.ARG.L); 
SPTCJ3  :*  OERIVCT.N.CT.NRl, O.ARG.L); 
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Sourc*  Listing 


lO-Osc-1984  16:37151 
lO-Dsc-1984  16:37:39 


VAX—] 

_0«AC 


J :  *  j  ♦  1  • 

end: 

end: 

WRITELNC  *00  TOU  WANT  PLOTING  YSS-1  *, 
'PRINT-2  COHPARE-3  N0-4')l 
REA0CANS2); 

IF  ANS2-  1  THEN  BEGIN 
INITTC960): 

VUINDOCO. 0.2 56.0. 0.0. 256.0)1 
hovea<xco3.tco3)  : 

FOR  I :»0  TO  M-l  00 
DRAWACXCI3.YCI3): 

FOR  11-0  TO  Jl-1  00  BEGIN 

M0VEA(SPX1CI3-1.5,SPT1 Cl  3-1.5); 
ANCH0C111): 

end: 

MOVEA{SPXC03,SPTC03): 

FOR  IS*.  TO  J-l  00 

0ASHA(SPXCI3.SPTCI3.2): 

F1NITTC0.0): 

ENO 

ELSE  if  ANS2-2  THEN  BEGIN 
FOR  l:«0  TO  Jl-1  DO 

WRITELN«'SPX1-',SPX1CI3. 

*  SPn»',sPrici3): 

ENO 

ELSE  IF  ANS2-3  THEN  BEGIN 
WRITELNE'TOL*'); 

RE  AO CT OLE); 

splcoef; 

aluhp; 

ak :» i ;  RA:»XCH-13-XC03J 
WHILE  ENC  AK  3  >0  00  BEGIN 
TEHPXI-BEGCAKJ :  POSINX; 

starcak3:><<tehpt-xco3)-Ra/2)/ra: 
TEHPx:»ENCAK3:  POSINX: 
TERCAX3:»<(TEHPT-XC03)-RA/2)/RA; 

tehpx:»haxctcax3:  posinx; 
PEAKCAX3:*CCTEHPT— XI03 )-RA/2)/RA; 
AREACAK3  :><AREAC AK3)/CRA«*2 >: 
WRITELN('BEGIN-',STARCJK3,'  ENO*', 
TERC AK3. '  TOTAL-'. RA. 

*  AR6A«',AREACAK3, 

'  peak-'.peakcak3); 
ak:-ak*i; 

ENO; 

eno: 

WRITELNC'DO  TOU  WANT  RUN  AGAIN  YES-1'. 
'  NO-2'); 

reac(ans); 

IF  ANS-1  THEN  BEGIN 
WRITELNC'IOPT-') : 

REACCIOPT): 

eno; 

eno;  (*  while  •) 

ENO. 
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